Historical Contingency of Host-Microbiome Response to Perturbation

Microbiome Results

Overview

Setup

(hide)

Click on tabs to display additional information.

Libraries

# Load required libraries for statistical analysis and table creation
library(dplyr)
library(ggplot2)
library(gt)
library(broom)
library(car)
library(emmeans)
library(multcomp)
library(MASS)  # For negative binomial regression
library(rstatix)
library(coin)
library(phyloseq)
library(microViz)
library(tidyverse)

Plotting

# Define treatment order and color palette
treatment_order <- c(
  "A- T- P-",  # Control
  "A- T- P+",  # Parasite
  "A+ T- P-",  # Antibiotics
  "A+ T- P+",  # Antibiotics_Parasite
  "A- T+ P-",  # Temperature
  "A- T+ P+",  # Temperature_Parasite
  "A+ T+ P-",  # Antibiotics_Temperature
  "A+ T+ P+"   # Antibiotics_Temperature_Parasite
)

# Custom color palette matching treatment order
treatment_colors <- c(
  "#1B9E77",  # A- T- P- (Control)
  "#D95F02",  # A- T- P+ (Parasite)
  "#7570B3",  # A+ T- P- (Antibiotics)
  "#E7298A",  # A+ T- P+ (Antibiotics_Parasite)
  "#66A61E",  # A- T+ P- (Temperature)
  "#E6AB02",  # A- T+ P+ (Temperature_Parasite)
  "#A6761D",  # A+ T+ P- (Antibiotics_Temperature)
  "#666666"   # A+ T+ P+ (Antibiotics_Temperature_Parasite)
)

# Create named vector for color scale
treatment_color_scale <- setNames(treatment_colors, treatment_order)

Functions

# Function to extract sample data as dataframe from phyloseq object
samdatAsDataframe <- function(ps) {
  samdat <- phyloseq::sample_data(ps)
  df <- data.frame(samdat, check.names = FALSE, stringsAsFactors = FALSE)
  return(df)
}

# Function to rename variables in phyloseq object
ps_rename <- function(ps, ...) {
  ps <- microViz::ps_get(ps)
  df <- samdatAsDataframe(ps)
  df <- dplyr::rename(.data = df, ...)
  phyloseq::sample_data(ps) <- df
  return(ps)
}

# SourceFolder function
source(here::here("Code", "R", "Functions", "StartFunctions", "sourceFolder.R"))

# Import all helper functions found in `/Functions`
sourceFolder(here::here("Code", "R", "Functions", "StartFunctions"), T)

## 
## Attaching package: 'scales'

## The following object is masked from 'package:purrr':
## 
##     discard

## The following object is masked from 'package:readr':
## 
##     col_factor

## The following object is masked from 'package:coin':
## 
##     pvalue

## Loading required package: future

## 
## Attaching package: 'future'

## The following object is masked from 'package:survival':
## 
##     cluster

## 
## Attaching package: 'ggpubr'

## The following object is masked from 'package:microViz':
## 
##     stat_chull

## Loading required package: grid

## ========================================
## ComplexHeatmap version 2.18.0
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
## 
## If you use it in published research, please cite either one:
## - Gu, Z. Complex Heatmap Visualization. iMeta 2022.
## - Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
##     genomic data. Bioinformatics 2016.
## 
## 
## The new InteractiveComplexHeatmap package can directly export static 
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
## 
## This message can be suppressed by:
##   suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================

## Loading required package: data.table

## 
## Attaching package: 'data.table'

## The following objects are masked from 'package:lubridate':
## 
##     hour, isoweek, mday, minute, month, quarter, second, wday, week,
##     yday, year

## The following object is masked from 'package:purrr':
## 
##     transpose

## The following objects are masked from 'package:dplyr':
## 
##     between, first, last

## Loading required package: doParallel

## Loading required package: foreach

## 
## Attaching package: 'foreach'

## The following objects are masked from 'package:purrr':
## 
##     accumulate, when

## Loading required package: iterators

## Loading required package: GUniFrac

## Registered S3 method overwritten by 'rmutil':
##   method         from
##   print.response httr

## Loading required package: magrittr

## 
## Attaching package: 'magrittr'

## The following object is masked from 'package:purrr':
## 
##     set_names

## The following object is masked from 'package:tidyr':
## 
##     extract

## Loading required package: rlang

## 
## Attaching package: 'rlang'

## The following object is masked from 'package:magrittr':
## 
##     set_names

## The following object is masked from 'package:data.table':
## 
##     :=

## The following objects are masked from 'package:purrr':
## 
##     %@%, flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl,
##     flatten_raw, invoke, splice

## Loading required package: vegan

## Loading required package: permute

## Loading required package: lattice

## This is vegan 2.6-4

## 
## microbiome R package (microbiome.github.com)
##     
## 
## 
##  Copyright (C) 2011-2022 Leo Lahti, 
##     Sudarshan Shetty et al. <microbiome.github.io>

## 
## Attaching package: 'microbiome'

## The following object is masked from 'package:vegan':
## 
##     diversity

## The following object is masked from 'package:scales':
## 
##     alpha

## The following object is masked from 'package:ggplot2':
## 
##     alpha

## The following object is masked from 'package:base':
## 
##     transform

## Loading required package: ape

## 
## Attaching package: 'ape'

## The following object is masked from 'package:ggpubr':
## 
##     rotate

## The following object is masked from 'package:dplyr':
## 
##     where

## Loading required package: nlme

## 
## Attaching package: 'nlme'

## The following object is masked from 'package:dplyr':
## 
##     collapse

## Loading required package: Matrix

## 
## Attaching package: 'Matrix'

## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

## 
## Attaching package: 'lme4'

## The following object is masked from 'package:nlme':
## 
##     lmList

## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## 3 files sourced from folder "/Users/michaelsieler/Dropbox/Mac (2)/Documents/Sharpton_Lab/Projects_Repository/Rules_of_Life/major-experiment-2023/Code/R/Functions/StartFunctions"

sourceFolder(here::here("Code", "R", "Functions", "HelperFunctions"), T)

## 9 files sourced from folder "/Users/michaelsieler/Dropbox/Mac (2)/Documents/Sharpton_Lab/Projects_Repository/Rules_of_Life/major-experiment-2023/Code/R/Functions/HelperFunctions"

# sourceFolder(here::here("Code", "R", "Functions", "AnalysisScripts"), T)

# Source the calculate_dynamic_metrics.R script to load required functions
source(here::here("Code", "R", "Functions", "AnalysisScripts", "calculate_dynamic_metrics.R"))

Import Data

ps.tmp <- readRDS("/Users/michaelsieler/Dropbox/Mac (2)/Documents/Sharpton_Lab/Projects_Repository/Rules_of_Life/major-experiment-2023/Data/Robjects/pseq_uncleaned_05052025.rds")


ps.cleaned <-
    ps.tmp %>%
        ## Update Metadata
        ps_rename(Time = Timepoint) %>%
        microViz::ps_mutate(
            Treatment = case_when(
                Antibiotics == 0 & Temperature == 0 & Pathogen == 0 ~ "A- T- P-",
                Antibiotics == 0 & Temperature == 0 & Pathogen == 1 ~ "A- T- P+",
                Antibiotics == 1 & Temperature == 0 & Pathogen == 0 ~ "A+ T- P-",
                Antibiotics == 1 & Temperature == 0 & Pathogen == 1 ~ "A+ T- P+",
                Antibiotics == 0 & Temperature == 1 & Pathogen == 0 ~ "A- T+ P-",
                Antibiotics == 0 & Temperature == 1 & Pathogen == 1 ~ "A- T+ P+",
                Antibiotics == 1 & Temperature == 1 & Pathogen == 0 ~ "A+ T+ P-",
                Antibiotics == 1 & Temperature == 1 & Pathogen == 1 ~ "A+ T+ P+",
                TRUE ~ "Unknown"
            ), .after = "Pathogen"
        ) %>%
        microViz::ps_mutate(Sample = fecal.sample.number, .before = 1) %>%
        microViz::ps_mutate(Sample = gsub("^f", "", Sample)) %>%
        microViz::ps_filter(Treatment != "Unknown") %>%
        microViz::ps_mutate(
            History = case_when(
                Antibiotics + Temperature == 0 ~ 0,
                Antibiotics + Temperature == 1 ~ 1,
                Antibiotics + Temperature == 2 ~ 2,
            ), .after = "Treatment"
        ) %>%
        
        ## Additional metadata updates, factorizing metadata
        microViz::ps_mutate(
        # Create treatment code
            treatment_code = case_when(
              Antibiotics == 0 & Temperature == 0 & Pathogen == 0 ~ "Aneg_Tneg_Pneg",
              Antibiotics == 0 & Temperature == 0 & Pathogen == 1 ~ "Aneg_Tneg_Ppos",
              Antibiotics == 1 & Temperature == 0 & Pathogen == 0 ~ "Apos_Tneg_Pneg",
              Antibiotics == 1 & Temperature == 0 & Pathogen == 1 ~ "Apos_Tneg_Ppos",
              Antibiotics == 0 & Temperature == 1 & Pathogen == 0 ~ "Aneg_Tpos_Pneg",
              Antibiotics == 0 & Temperature == 1 & Pathogen == 1 ~ "Aneg_Tpos_Ppos",
              Antibiotics == 1 & Temperature == 1 & Pathogen == 0 ~ "Apos_Tpos_Pneg",
              Antibiotics == 1 & Temperature == 1 & Pathogen == 1 ~ "Apos_Tpos_Ppos"
            ),
            # Create treatment group factor
            treatment_group = case_when(
              Antibiotics == 0 & Temperature == 0 & Pathogen == 1 ~ "Parasite",
              Antibiotics == 1 & Temperature == 0 & Pathogen == 0 ~ "Antibiotics",
              Antibiotics == 1 & Temperature == 0 & Pathogen == 1 ~ "Antibiotics_Parasite",
              Antibiotics == 0 & Temperature == 1 & Pathogen == 0 ~ "Temperature",
              Antibiotics == 0 & Temperature == 1 & Pathogen == 1 ~ "Temperature_Parasite",
              Antibiotics == 1 & Temperature == 1 & Pathogen == 0 ~ "Antibiotics_Temperature",
              Antibiotics == 1 & Temperature == 1 & Pathogen == 1 ~ "Antibiotics_Temperature_Parasite",
              TRUE ~ "Control"
            ),
            # Convert to factor with appropriate levels
            treatment_group = factor(treatment_group, 
                                   levels = c("Control", "Parasite", 
                                              "Antibiotics", "Antibiotics_Parasite",
                                              "Temperature", "Temperature_Parasite",
                                            "Antibiotics_Temperature", "Antibiotics_Temperature_Parasite")
                                   ),
            treatment_code = factor(treatment_code, levels = treatment_order),
            # Create time point factor
            time_point = factor(Time, levels = c(0, 14, 18, 25, 29, 60)),
            # Create pathogen status factor
            pathogen_status = factor(ifelse(Pathogen == 1, "Exposed", "Unexposed"),
                                   levels = c("Unexposed", "Exposed")),
            # Create sex factor
            sex = factor(Sex, levels = c("M", "F"))
            )  %>%
    microViz::ps_mutate(Treatment = factor(Treatment, levels = treatment_order)) %>%
      microViz::ps_mutate(Exp_Type = case_when(
          Treatment %in% c("A- T- P-", "A- T- P+")  ~ "No prior stressor(s)",
          Treatment %in% c("A+ T- P-", "A+ T- P+")  ~ "Antibiotics",
          Treatment %in% c("A- T+ P-", "A- T+ P+") ~ "Temperature",
          Treatment %in% c("A+ T+ P-", "A+ T+ P+") ~ "Combined",
      )) %>%
      microViz::ps_mutate(Exp_Type = factor(Exp_Type, levels = c("No prior stressor(s)", "Antibiotics", "Temperature", "Combined"))) %>%
  # Fix names for taxonomic ranks not identified
  microViz::tax_fix(suffix_rank = "current", anon_unique = T, unknown = NA) %>% 
  # Filter for any samples that contain more than 5000 reads
  microViz::ps_filter(sample_sums(.) > 5000) %>%
  # Any taxa not found in at least 3 samples are removed
  microViz::tax_filter(min_prevalence = 3, undetected = 0) %>%
  # Remove any unwanted reads
  microViz::tax_select(c("Mitochondria", "Chloroplast", "Eukaryota"), deselect = TRUE) %>%
  microViz::tax_select(c("Bacteria, Phylum"), deselect = TRUE) 


# ps.cleaned %>% microViz::samdat_tbl()

Variables

## Diversity Methods

diversity.method <- list()
diversity.method[["alpha"]] <- c("shannon", 
                                 "inverse_simpson", 
                                 "observed")
diversity.method[["beta"]] <- c("bray", "canberra")

## Stats/Plotting Variables

alpha.stats <- list() # save all stat results here
alpha.plots <- list() # save plots here

beta.stats <- list() # save all stat results here
beta.plots <- list() # save plots here

diffAbnd.stats <- list() # save all stat results here
diffAbnd.plots <- list() # save plots here

Calculate Scores

Alpha

## Alpha
ps.cleaned_alpha <-
    ps.cleaned %>% # phyloseq object we'll be using
  microViz::ps_calc_diversity( # calculates various alpha diversity metrics
    rank = "Genus", # What taxonomic rank do you want to calculate diversity metrics at
    index = "shannon", # What diversity metric to use
    varname = "Shannon__Genus", # Column name in sample data
    exp = F # exponentiate the result or not
  ) %>%
  microViz::ps_calc_diversity(
    rank = "Genus",
    index = "gini_simpson",
    varname = "Simpson__Genus" # Column name in sample data
  ) %>%
  microViz::ps_calc_richness( # related to ps_calc_diversity, this calculates richness or "observed" values
    rank = "Genus",
    varname = "Richness__Genus"
  ) %>%
  # ps_calc_diversity.phy(
  #   varname = "Phylogenetic__Genus"
  # ) %>% # Note: This is a helper function I made. You need to adjust function manually to change taxon rank (See: microViz_Helper.R)
  # # ps_calc_diversity.phy() helper function does not properly name the column. Modify the following function as needed 
  # ps_rename("Phylogenetic__Genus" = "phylogenetic_Genus") %>% # Helper function (See: microViz_Helper.R)
  ps_mutate(across(contains("__Genus"), norm_scores, .names = "{.col}_norm")) 

## 
##     lambda     W Shapiro.p.value    A Anderson.p.value
## 420  0.475 0.987       3.787e-05 2.49          2.7e-06
## 
## if (lambda >  0){TRANS = x ^ lambda} 
## if (lambda == 0){TRANS = log(x)} 
## if (lambda <  0){TRANS = -1 * x ^ lambda} 
## 
## 
##     lambda      W Shapiro.p.value     A Anderson.p.value
## 482  2.025 0.9761       2.845e-08 2.862         3.33e-07
## 
## if (lambda >  0){TRANS = x ^ lambda} 
## if (lambda == 0){TRANS = log(x)} 
## if (lambda <  0){TRANS = -1 * x ^ lambda} 
## 
## 
##     lambda      W Shapiro.p.value    A Anderson.p.value
## 393   -0.2 0.9883       0.0001119 1.69        0.0002453
## 
## if (lambda >  0){TRANS = x ^ lambda} 
## if (lambda == 0){TRANS = log(x)} 
## if (lambda <  0){TRANS = -1 * x ^ lambda}

ps.cleaned_alpha %>%
    microViz::samdat_tbl()

## # A tibble: 594 × 37
##    .sample_name Sample Tank.ID  Time Antibiotics Temperature Pathogen Treatment
##    <chr>        <chr>    <dbl> <dbl>       <dbl>       <dbl>    <dbl> <fct>    
##  1 f1           1            1     0           0           0        0 A- T- P- 
##  2 f2           2            1     0           0           0        0 A- T- P- 
##  3 f3           3            1     0           0           0        0 A- T- P- 
##  4 f31          31           2     0           0           0        0 A- T- P- 
##  5 f32          32           2     0           0           0        0 A- T- P- 
##  6 f33          33           2     0           0           0        0 A- T- P- 
##  7 f61          61           3     0           0           0        0 A- T- P- 
##  8 f62          62           3     0           0           0        0 A- T- P- 
##  9 f63          63           3     0           0           0        0 A- T- P- 
## 10 f91          91           4     0           0           0        1 A- T- P+ 
## # ℹ 584 more rows
## # ℹ 29 more variables: History <dbl>, Sex <chr>, Weight.g <dbl>,
## #   Weight.mg <dbl>, Length.mm <dbl>, Body.Size <chr>, Mortality <dbl>,
## #   Total.Worm.Count <dbl>, Flagged <dbl>, notes <chr>,
## #   fecal.sample.number <chr>, gut.sample.number <chr>, adultworm.count <dbl>,
## #   larvalworm.count <dbl>, notes.dissection <chr>, PLATE <chr>, WELL <chr>,
## #   treatment_code <fct>, treatment_group <fct>, time_point <fct>, …

test.data.alpha <-
    ps.cleaned_alpha %>% 
    microViz::samdat_tbl() %>%
    dplyr::relocate(c(".sample_name"), .before = "Antibiotics") %>%
    tidyr::pivot_longer(cols = contains("__Genus"), names_to = "Alpha.Metric", values_to = "Alpha.Score") %>%
    dplyr::mutate(Data.Type = case_when(
        Alpha.Metric %in% c("Shannon__Genus_norm", "Simpson__Genus_norm", "Richness__Genus_norm") ~ "Normalized",
        TRUE ~ "Raw"
    )) %>%
    dplyr::relocate(c("Alpha.Metric", "Alpha.Score", "Data.Type"), .after = "Tank.ID") %>%
    dplyr::relocate(c("Sample", "Treatment", "Time", "Tank.ID"), .before = 1)
    # tidyr::pivot_longer(cols = contains("_norm"), names_to = "Alpha.Metric", values_to = "Alpha.Score") %>%
    # dplyr::filter(Alpha.Metric %in% c("Simpson__Genus_norm", "Simpson__Genus_norm")) %>%
    # tidyr::pivot_longer(cols = contains("__Genus"), names_to = "Raw.Metric", values_to = "Raw.Score") 

Beta

Dist Matrices

#### Beta diversity matrices

##### All

beta.dist.mat <- list()

beta.dist.mat[["All"]] <- # saves the results of this loop to a list
  furrr::future_map(diversity.method[["beta"]], function(x){ # Future_map runs the function in parallel
    ps.cleaned_alpha %>% # is the phyloseq object we are looping over
      microViz::tax_transform(trans = "identity", # tax transform transforms taxa counts
                    rank = ifelse(x == "gunifrac" || x == "wunifrac" || x == "unifrac", "unique", "Genus")) %>% # phylogenetic methods need "unique" ASV counts
      microViz::dist_calc(x, gunifrac_alpha = 0.5) # calculates the distance between samples. gunifrac_alpha = 0.5 weights abundance into its calculation
  }) %>%
  setNames(diversity.method[["beta"]]) # assigns names of the beta methods to the list

beta.dist.mat[["All_60DPE"]] <- # saves the results of this loop to a list
  furrr::future_map(diversity.method[["beta"]], function(x){ # Future_map runs the function in parallel
    ps.cleaned_alpha %>% # is the phyloseq object we are looping over
          microViz::ps_filter(Time == 60) %>%
      microViz::tax_transform(trans = "identity", # tax transform transforms taxa counts
                    rank = ifelse(x == "gunifrac" || x == "wunifrac" || x == "unifrac", "unique", "Genus")) %>% # phylogenetic methods need "unique" ASV counts
      microViz::dist_calc(x, gunifrac_alpha = 0.5) # calculates the distance between samples. gunifrac_alpha = 0.5 weights abundance into its calculation
  }) %>%
  setNames(diversity.method[["beta"]]) # assigns names of the beta methods to the list

##### Parasite exposure only

beta.dist.mat[["ParasiteExp"]] <- # saves the results of this loop to a list
  furrr::future_map(diversity.method[["beta"]], function(x){ # Future_map runs the function in parallel
    ps.cleaned_alpha %>% # is the phyloseq object we are looping over
          microViz::ps_filter(Treatment  %in% c(
          "A- T- P-",
          "A- T- P+"#,
          # "A+ T- P-",
          # "A+ T- P+",
          # "A- T+ P-",
          # "A- T+ P+",
          # "A+ T+ P-",
          # "A+ T+ P+"
          )) %>%
      microViz::tax_transform(trans = "identity", # tax transform transforms taxa counts
                    rank = ifelse(x == "gunifrac" || x == "wunifrac" || x == "unifrac", "unique", "Genus")) %>% # phylogenetic methods need "unique" ASV counts
      microViz::dist_calc(x, gunifrac_alpha = 0.5) # calculates the distance between samples. gunifrac_alpha = 0.5 weights abundance into its calculation
  }) %>%
  setNames(diversity.method[["beta"]]) # assigns names of the beta methods to the list

beta.dist.mat[["ParasiteExp_60DPE"]] <- # saves the results of this loop to a list
  furrr::future_map(diversity.method[["beta"]], function(x){ # Future_map runs the function in parallel
    ps.cleaned_alpha %>% # is the phyloseq object we are looping over
          microViz::ps_filter(Treatment  %in% c(
          "A- T- P-",
          "A- T- P+"#,
          # "A+ T- P-",
          # "A+ T- P+",
          # "A- T+ P-",
          # "A- T+ P+",
          # "A+ T+ P-",
          # "A+ T+ P+"
          )) %>%
          microViz::ps_filter(Time == 60) %>%
      microViz::tax_transform(trans = "identity", # tax transform transforms taxa counts
                    rank = ifelse(x == "gunifrac" || x == "wunifrac" || x == "unifrac", "unique", "Genus")) %>% # phylogenetic methods need "unique" ASV counts
      microViz::dist_calc(x, gunifrac_alpha = 0.5) # calculates the distance between samples. gunifrac_alpha = 0.5 weights abundance into its calculation
  }) %>%
  setNames(diversity.method[["beta"]]) # assigns names of the beta methods to the list

##### Prior stressors and parasite exposure

beta.dist.mat[["PriorStressParaExp"]] <- # saves the results of this loop to a list
  furrr::future_map(diversity.method[["beta"]], function(x){ # Future_map runs the function in parallel
    ps.cleaned_alpha %>% # is the phyloseq object we are looping over
          microViz::ps_filter(Treatment  %in% c(
          # "A- T- P-",
          "A- T- P+",
          # "A+ T- P-",
          "A+ T- P+",
          # "A- T+ P-",
          "A- T+ P+",
          # "A+ T+ P-",
          "A+ T+ P+"
          )) %>%
      microViz::tax_transform(trans = "identity", # tax transform transforms taxa counts
                    rank = ifelse(x == "gunifrac" || x == "wunifrac" || x == "unifrac", "unique", "Genus")) %>% # phylogenetic methods need "unique" ASV counts
      microViz::dist_calc(x, gunifrac_alpha = 0.5) # calculates the distance between samples. gunifrac_alpha = 0.5 weights abundance into its calculation
  }) %>%
  setNames(diversity.method[["beta"]]) # assigns names of the beta methods to the list

beta.dist.mat[["PriorStressParaExp_60DPE"]] <- # saves the results of this loop to a list
  furrr::future_map(diversity.method[["beta"]], function(x){ # Future_map runs the function in parallel
    ps.cleaned_alpha %>% # is the phyloseq object we are looping over
          microViz::ps_filter(Treatment  %in% c(
          # "A- T- P-",
          "A- T- P+",
          # "A+ T- P-",
          "A+ T- P+",
          # "A- T+ P-",
          "A- T+ P+",
          # "A+ T+ P-",
          "A+ T+ P+"
          )) %>%
          microViz::ps_filter(Time == 60) %>%
      microViz::tax_transform(trans = "identity", # tax transform transforms taxa counts
                    rank = ifelse(x == "gunifrac" || x == "wunifrac" || x == "unifrac", "unique", "Genus")) %>% # phylogenetic methods need "unique" ASV counts
      microViz::dist_calc(x, gunifrac_alpha = 0.5) # calculates the distance between samples. gunifrac_alpha = 0.5 weights abundance into its calculation
  }) %>%
  setNames(diversity.method[["beta"]]) # assigns names of the beta methods to the list

##### Prior stressors and NO parasite exposure

beta.dist.mat[["PriorStressNoParaExp"]] <- # saves the results of this loop to a list
  furrr::future_map(diversity.method[["beta"]], function(x){ # Future_map runs the function in parallel
    ps.cleaned_alpha %>% # is the phyloseq object we are looping over
          microViz::ps_filter(Treatment  %in% c(
          "A- T- P-",
          # "A- T- P+",
          "A+ T- P-",
          # "A+ T- P+",
          "A- T+ P-",
          # "A- T+ P+",
          "A+ T+ P-"#,
          # "A+ T+ P+"
          )) %>%
      microViz::tax_transform(trans = "identity", # tax transform transforms taxa counts
                    rank = ifelse(x == "gunifrac" || x == "wunifrac" || x == "unifrac", "unique", "Genus")) %>% # phylogenetic methods need "unique" ASV counts
      microViz::dist_calc(x, gunifrac_alpha = 0.5) # calculates the distance between samples. gunifrac_alpha = 0.5 weights abundance into its calculation
  }) %>%
  setNames(diversity.method[["beta"]]) # assigns names of the beta methods to the list

beta.dist.mat[["PriorStressNoParaExp_60DPE"]] <- # saves the results of this loop to a list
  furrr::future_map(diversity.method[["beta"]], function(x){ # Future_map runs the function in parallel
    ps.cleaned_alpha %>% # is the phyloseq object we are looping over
          microViz::ps_filter(Treatment  %in% c(
          "A- T- P-",
          # "A- T- P+",
          "A+ T- P-",
          # "A+ T- P+",
          "A- T+ P-",
          # "A- T+ P+",
          "A+ T+ P-"#,
          # "A+ T+ P+"
          )) %>%
          microViz::ps_filter(Time == 60) %>%
      microViz::tax_transform(trans = "identity", # tax transform transforms taxa counts
                    rank = ifelse(x == "gunifrac" || x == "wunifrac" || x == "unifrac", "unique", "Genus")) %>% # phylogenetic methods need "unique" ASV counts
      microViz::dist_calc(x, gunifrac_alpha = 0.5) # calculates the distance between samples. gunifrac_alpha = 0.5 weights abundance into its calculation
  }) %>%
  setNames(diversity.method[["beta"]]) # assigns names of the beta methods to the list

Personalization

#### Beta Personalization Scores

data.beta.list <- list()

# Save dataframe
data.beta.list[["All"]] <-
  
  ps.cleaned_alpha %>% 
  
  # Calculate Distance matrix
  microViz::dist_calc("bray") %>% 
  
  # Get distance matric
  microViz::dist_get() %>%
  # Pat Schloss methods to view beta diversity matrix data
  # Source: https://youtu.be/Gm-kg3TuML4?si=nYLLAxm7uv8aYG27&t=250
  
  # Convert into a matrix object
  as.matrix() %>%
  
  # Convert matrix into a tibble
  tibble::as_tibble(rownames = "Sample_a") %>%
  
  # Convert dataframe into "long" format
  tidyr::pivot_longer(-Sample_a, names_to = "Sample_b", values_to = "Dist") %>%
  
  # Remove any rows where names equal eachother
  dplyr::filter(Sample_a != Sample_b) %>%
  
  # Merge sample dataframe
  dplyr::right_join((ps.cleaned_alpha %>% 
                       microViz::samdat_tbl() %>% ungroup() %>%  arrange(Time, Treatment, Tank.ID) %>% 
                       dplyr::select(!contains("_Genus")) %>%
                       dplyr::select(fecal.sample.number, Time, Treatment, Tank.ID)
  ), 
  by = c("Sample_a" = "fecal.sample.number")) %>% 
  dplyr::right_join((ps.cleaned_alpha %>% 
                       microViz::samdat_tbl() %>% 
                       dplyr::ungroup() %>% 
                       dplyr::arrange(Time, Treatment, Tank.ID) %>% 
                       dplyr::select(!contains("_Genus")) %>%
                       dplyr::select(fecal.sample.number, Time, Treatment, Tank.ID)
  ), 
  by = c("Sample_b" = "fecal.sample.number"), suffix = c("_a", "_b")) %>% # add suffix to differentiate between sample metadata
  # Move columns around
  dplyr::select(order(colnames(.))) %>%
  dplyr::relocate(dplyr::any_of(c("Sample_a", "Sample_b", "Dist")), .before = 1) %>%
  dplyr::relocate(dplyr::contains(c("Tank")), .after = dplyr::last_col()) %>%
  dplyr::ungroup() %>%
  
  # Calculate the median distance of samples to all other samples (Measure of "personalization")
  dplyr::group_by(Sample_a, Treatment_a, Time_a, Tank.ID_a) %>%
  dplyr::summarize(Beta.Score = median(Dist), #*100,  
                   Beta.Score.Mean = mean(Dist), #*100, 
                   count = dplyr::n()) %>% 
  dplyr::ungroup() %>%
  
  # Rename Columns
  dplyr::rename(fecal.sample.number = Sample_a,
                Treatment = Treatment_a,
                
                Time = Time_a,
                Tank.ID = Tank.ID_a) %>%
  dplyr::right_join(ps.cleaned_alpha %>%
                      microViz::samdat_tbl() %>% 
                      dplyr::select(!contains("_Genus")),
                    by = c("fecal.sample.number", "Time", "Treatment","Tank.ID")) %>%
  dplyr::ungroup() %>%
  
  # Normalize diversity scores
  dplyr::mutate(dplyr::across(dplyr::contains("Beta.Score"), norm_scores, .names = "{.col}_norm"), .after = "Beta.Score") %>%
  
  # Create a column to save what type of beta metric was calculated
  dplyr::mutate(Beta.Metric = "Bray-Curtis", .before = "Beta.Score") %>%
  dplyr::mutate(Sample = gsub("^f", "", fecal.sample.number)) %>%
    dplyr::relocate(Sample, .before = 1) %>%
    dplyr::relocate(fecal.sample.number, .after = 36)

## `summarise()` has grouped output by 'Sample_a', 'Treatment_a', 'Time_a'. You
## can override using the `.groups` argument.

## 
##     lambda      W Shapiro.p.value    A Anderson.p.value
## 291  -2.75 0.9882       0.0001004 1.05         0.009232
## 
## if (lambda >  0){TRANS = x ^ lambda} 
## if (lambda == 0){TRANS = log(x)} 
## if (lambda <  0){TRANS = -1 * x ^ lambda} 
## 
## 
##     lambda      W Shapiro.p.value     A Anderson.p.value
## 235  -4.15 0.9883       0.0001067 1.004          0.01196
## 
## if (lambda >  0){TRANS = x ^ lambda} 
## if (lambda == 0){TRANS = log(x)} 
## if (lambda <  0){TRANS = -1 * x ^ lambda}

data.beta.list[["All"]] <- data.beta.list[["All"]] %>%
dplyr::full_join(., 
                   ps.cleaned_alpha %>% 
                     # tax_agg("Genus") %>%
                     microViz::dist_calc("canberra") %>% 
                     
                     # Get distance matric
    microViz::dist_get() %>%
    # Pat Schloss methods to view beta diversity matrix data
    # Source: https://youtu.be/Gm-kg3TuML4?si=nYLLAxm7uv8aYG27&t=250
    
    # Convert into a matrix object
    as.matrix() %>%
    
    # Convert matrix into a tibble
    tibble::as_tibble(rownames = "Sample_a") %>%
    
    # Convert dataframe into "long" format
    tidyr::pivot_longer(-Sample_a, names_to = "Sample_b", values_to = "Dist") %>%
    
    # Remove any rows where names equal eachother
    dplyr::filter(Sample_a != Sample_b) %>%
    
    # Merge sample dataframe
    dplyr::right_join((ps.cleaned_alpha %>% 
                       microViz::samdat_tbl() %>% ungroup() %>%  arrange(Time, Treatment, Tank.ID) %>% 
                       dplyr::select(!contains("_Genus")) %>%
                       dplyr::select(fecal.sample.number, Time, Treatment, Tank.ID)
    ), 
    by = c("Sample_a" = "fecal.sample.number")) %>% 
    dplyr::right_join((ps.cleaned_alpha %>% 
                       microViz::samdat_tbl() %>% 
                       dplyr::ungroup() %>% 
                       dplyr::arrange(Time, Treatment, Tank.ID) %>% 
                       dplyr::select(!contains("_Genus")) %>%
                       dplyr::select(fecal.sample.number, Time, Treatment, Tank.ID)
    ), 
    by = c("Sample_b" = "fecal.sample.number"), suffix = c("_a", "_b")) %>% # add suffix to differentiate between sample metadata
    # Move columns around
    dplyr::select(order(colnames(.))) %>%
    dplyr::relocate(dplyr::any_of(c("Sample_a", "Sample_b", "Dist")), .before = 1) %>%
    dplyr::relocate(dplyr::contains(c("Tank")), .after = dplyr::last_col()) %>%
    dplyr::ungroup() %>%
    
    # Calculate the median distance of samples to all other samples (Measure of "personalization")
    dplyr::group_by(Sample_a, Treatment_a, Time_a, Tank.ID_a) %>%
    dplyr::summarize(Beta.Score = median(Dist), #*100,  
                   Beta.Score.Mean = mean(Dist), #*100, 
                   count = dplyr::n()) %>% 
    dplyr::ungroup() %>%
    
    # Rename Columns
    dplyr::rename(fecal.sample.number = Sample_a,
                Treatment = Treatment_a,
                
                Time = Time_a,
                Tank.ID = Tank.ID_a) %>%
    dplyr::right_join(ps.cleaned_alpha %>%
                      microViz::samdat_tbl() %>% 
                      dplyr::select(!contains("_Genus")),
                    by = c("fecal.sample.number", "Time", "Treatment","Tank.ID")) %>%
    dplyr::ungroup() %>%
                     
    # Normalize diversity scores
    dplyr::mutate(dplyr::across(contains("Beta.Score"), norm_scores, .names = "{.col}_norm"), .after = "Beta.Score") %>%
    
    # Create a column to save what type of beta metric was calculated
    dplyr::mutate(Beta.Metric = "Canberra", .before = "Beta.Score")
  )

## `summarise()` has grouped output by 'Sample_a', 'Treatment_a', 'Time_a'. You
## can override using the `.groups` argument.

## 
##    lambda      W Shapiro.p.value      A Anderson.p.value
## 76 -8.125 0.9925        0.004258 0.6477          0.09061
## 
## if (lambda >  0){TRANS = x ^ lambda} 
## if (lambda == 0){TRANS = log(x)} 
## if (lambda <  0){TRANS = -1 * x ^ lambda} 
## 
## 
##    lambda      W Shapiro.p.value      A Anderson.p.value
## 51  -8.75 0.9948         0.04029 0.4685           0.2483
## 
## if (lambda >  0){TRANS = x ^ lambda} 
## if (lambda == 0){TRANS = log(x)} 
## if (lambda <  0){TRANS = -1 * x ^ lambda}

## Joining with `by = join_by(Sample, Treatment, Time, Tank.ID, Beta.Metric,
## Beta.Score, Beta.Score_norm, Beta.Score.Mean_norm, Beta.Score.Mean, count,
## .sample_name, Antibiotics, Temperature, Pathogen, History, Sex, Weight.g,
## Weight.mg, Length.mm, Body.Size, Mortality, Total.Worm.Count, Flagged, notes,
## gut.sample.number, adultworm.count, larvalworm.count, notes.dissection, PLATE,
## WELL, treatment_code, treatment_group, time_point, pathogen_status, sex,
## fecal.sample.number, Exp_Type)`

# Prepare beta diversity data
test.data.beta <-
    data.beta.list[["All"]] %>%
        tidyr::pivot_longer(cols = c("Beta.Score", "Beta.Score_norm"), names_to = "Beta.Metric.Type", values_to = "Beta.Score") %>%
        dplyr::relocate(c("Beta.Metric.Type", "Beta.Score"), .after = "Beta.Metric") %>%
        dplyr::mutate(Data.Type = case_when(
            Beta.Metric.Type %in% c("   Beta.Score", "Beta.Score_norm") ~ "Normalized",
            TRUE ~ "Raw"
        ), .after = "Beta.Score") %>%
        # Remove the temporary Beta.Metric.Type column
        dplyr::select(-c(Beta.Metric.Type, count)) %>%
        dplyr::mutate(Beta.Metric = ifelse(Data.Type == "Normalized", paste0(Beta.Metric, "_norm"), Beta.Metric)) %>%
        dplyr::arrange(Sample, Time, Tank.ID, Data.Type)

Calculate SADMMs

# Calculate all metrics at once using the combined function
test.alpha.all_metrics <- calculate_all_dynamic_metrics(
  data = test.data.alpha,
  metric_col = "Alpha.Metric",
  score_col = "Alpha.Score",
  ref_group = NULL,  # Using all groups as reference
  time_col = "Time",
  group_col = "Treatment",
  use_mean = TRUE
)

## Joining with `by = join_by(Sample, Time, Treatment, Tank.ID, Alpha.Metric,
## Data.Type, Displacement)`
## Joining with `by = join_by(Sample, Time, Treatment, Tank.ID, Alpha.Metric,
## Data.Type)`

# Calculate all metrics at once using the combined function
test.beta.all_metrics <- calculate_all_dynamic_metrics(
  data = test.data.beta,
  metric_col = "Beta.Metric",
  score_col = "Beta.Score",
  ref_group = NULL,
  time_col = "Time",
  group_col = "Treatment",
  use_mean = TRUE
)

## Joining with `by = join_by(Sample, Time, Treatment, Tank.ID, Beta.Metric,
## Data.Type, Displacement)`
## Joining with `by = join_by(Sample, Time, Treatment, Tank.ID, Beta.Metric,
## Data.Type)`

# Merge Displacement data with Alpha.Metric and Alpha.Score
data.alpha.merged <-
    test.alpha.all_metrics$volatility %>%
        # dplyr::filter()
        # First pivot the SADMM metrics
        # dplyr::group_by(Data.Type) %>%
        tidyr::pivot_longer(
            cols = c("Alpha.Score", "Ref.Score", "Displacement", "Velocity", "Acceleration", "Volatility"),
            names_to = "METRIC",
            values_to = "SCORE"
        ) %>%
        dplyr::mutate(METRIC = factor(METRIC, levels = c("Alpha.Score", "Ref.Score", "Displacement", "Velocity", "Acceleration", "Volatility")),
                      Data.Type = factor(Data.Type, levels = c("Raw", "Normalized")))



# Merge Displacement data with Alpha.Metric and Alpha.Score
data.beta.merged <-
    test.beta.all_metrics$volatility %>%
        # dplyr::filter()
        # First pivot the SADMM metrics
        # dplyr::group_by(Data.Type) %>%
        tidyr::pivot_longer(
            cols = c("Beta.Score", "Ref.Score", "Displacement", "Velocity", "Acceleration", "Volatility"),
            names_to = "METRIC",
            values_to = "SCORE"
        ) %>%
        dplyr::mutate(METRIC = factor(METRIC, levels = c("Beta.Score", "Ref.Score", "Displacement", "Velocity", "Acceleration", "Volatility")),
                      Data.Type = factor(Data.Type, levels = c("Raw", "Normalized")))

00) All

Alpha

All

Displacement

Plot

Code

# Set the subsection for this analysis
tmp.resSubSection <- "All"

# GLM Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]] <- 
  test.alpha.all_metrics$volatility %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models(formula_str = "Displacement ~ Treatment", 
                   alpha_metric_col = "Alpha.Metric", 
                   alpha_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["ANOVA"]] <-
  run_glm_anova(alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]])

# Tukey Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["Tukey"]] <-
  test.alpha.all_metrics$volatility %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm(., "Displacement", "Alpha.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

Tables

GLM Results
glm(Displacement ~ Treatment); All treatments
term	estimate	std.error	statistic	p.value	p.adj.sig
Shannon
(Intercept)	0.078	0.019	3.995	<0.001	****
A- T- P+	−0.045	0.027	−1.627	0.104	ns
A+ T- P-	−0.051	0.028	−1.815	0.070	ns
A+ T- P+	−0.071	0.028	−2.584	0.010	**
A- T+ P-	−0.004	0.028	−0.151	≥0.25	ns
A- T+ P+	−0.092	0.028	−3.302	0.001	**
A+ T+ P-	−0.113	0.030	−3.799	<0.001	***
A+ T+ P+	−0.061	0.029	−2.149	0.032	*
Simpson
(Intercept)	0.042	0.028	1.527	0.127	ns
A- T- P+	−0.049	0.039	−1.238	0.216	ns
A+ T- P-	−0.034	0.040	−0.841	≥0.25	ns
A+ T- P+	−0.085	0.039	−2.165	0.031	*
A- T+ P-	0.035	0.041	0.874	≥0.25	ns
A- T+ P+	−0.154	0.040	−3.867	<0.001	***
A+ T+ P-	−0.181	0.042	−4.270	<0.001	****
A+ T+ P+	−0.091	0.041	−2.232	0.026	*
Richness
(Intercept)	0.221	0.024	9.249	<0.001	****
A- T- P+	0.009	0.034	0.265	≥0.25	ns
A+ T- P-	−0.049	0.035	−1.430	0.153	ns
A+ T- P+	−0.045	0.034	−1.325	0.186	ns
A- T+ P-	−0.016	0.035	−0.467	≥0.25	ns
A- T+ P+	0.016	0.034	0.456	≥0.25	ns
A+ T+ P-	0.029	0.037	0.783	≥0.25	ns
A+ T+ P+	−0.054	0.035	−1.533	0.126	ns

ANOVA of GLM
ANOVA(GLM(Displacement ~ Treatment), type = 2); All treatments
term	statistic	df	p.value	sig
Shannon
Treatment	25.086	7.000	<0.001	***
Simpson
Treatment	42.565	7.000	<0.001	****
Richness
Treatment	11.161	7.000	0.132	ns

Pairwise Tukey's HSD, p.adj: Dunnett
Tukey(Displacement ~ Treatment); All treatments
term	.y.	group1	group2	estimate	std.error	statistic	adj.p.value	Variable
Shannon
Treatment	Alpha.Score	A- T- P+	A- T- P-	−0.045	0.027	−1.627	≥0.25	Treatment
Treatment	Alpha.Score	A+ T- P-	A- T- P-	−0.051	0.028	−1.815	≥0.25	Treatment
Treatment	Alpha.Score	A+ T- P+	A- T- P-	−0.071	0.028	−2.584	0.161	Treatment
Treatment	Alpha.Score	A- T+ P-	A- T- P-	−0.004	0.028	−0.151	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A- T- P-	−0.092	0.028	−3.302	0.022	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T- P-	−0.113	0.030	−3.799	0.004	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T- P-	−0.061	0.029	−2.149	≥0.25	Treatment
Treatment	Alpha.Score	A+ T- P-	A- T- P+	−0.006	0.028	−0.220	≥0.25	Treatment
Treatment	Alpha.Score	A+ T- P+	A- T- P+	−0.027	0.028	−0.968	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P-	A- T- P+	0.040	0.028	1.421	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A- T- P+	−0.047	0.028	−1.695	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T- P+	−0.068	0.030	−2.298	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T- P+	−0.017	0.029	−0.582	≥0.25	Treatment
Treatment	Alpha.Score	A+ T- P+	A+ T- P-	−0.021	0.028	−0.730	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P-	A+ T- P-	0.047	0.029	1.608	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A+ T- P-	−0.041	0.028	−1.446	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A+ T- P-	−0.062	0.030	−2.056	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A+ T- P-	−0.010	0.029	−0.360	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P-	A+ T- P+	0.067	0.029	2.348	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A+ T- P+	−0.020	0.028	−0.729	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A+ T- P+	−0.042	0.030	−1.392	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A+ T- P+	0.010	0.029	0.353	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A- T+ P-	−0.088	0.029	−3.044	0.048	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T+ P-	−0.109	0.031	−3.549	0.009	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T+ P-	−0.057	0.029	−1.936	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T+ P+	−0.021	0.030	−0.706	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T+ P+	0.031	0.029	1.058	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A+ T+ P-	0.052	0.031	1.684	≥0.25	Treatment
Simpson
Treatment	Alpha.Score	A- T- P+	A- T- P-	−0.049	0.039	−1.238	≥0.25	Treatment
Treatment	Alpha.Score	A+ T- P-	A- T- P-	−0.034	0.040	−0.841	≥0.25	Treatment
Treatment	Alpha.Score	A+ T- P+	A- T- P-	−0.085	0.039	−2.165	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P-	A- T- P-	0.035	0.041	0.874	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A- T- P-	−0.154	0.040	−3.867	0.003	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T- P-	−0.181	0.042	−4.270	<0.001	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T- P-	−0.091	0.041	−2.232	≥0.25	Treatment
Treatment	Alpha.Score	A+ T- P-	A- T- P+	0.015	0.040	0.372	≥0.25	Treatment
Treatment	Alpha.Score	A+ T- P+	A- T- P+	−0.037	0.039	−0.935	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P-	A- T- P+	0.084	0.041	2.070	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A- T- P+	−0.105	0.040	−2.645	0.140	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T- P+	−0.133	0.042	−3.128	0.038	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T- P+	−0.042	0.041	−1.041	≥0.25	Treatment
Treatment	Alpha.Score	A+ T- P+	A+ T- P-	−0.052	0.040	−1.287	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P-	A+ T- P-	0.069	0.041	1.672	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A+ T- P-	−0.120	0.040	−2.963	0.060	Treatment
Treatment	Alpha.Score	A+ T+ P-	A+ T- P-	−0.148	0.043	−3.421	0.015	Treatment
Treatment	Alpha.Score	A+ T+ P+	A+ T- P-	−0.057	0.041	−1.381	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P-	A+ T- P+	0.121	0.041	2.962	0.061	Treatment
Treatment	Alpha.Score	A- T+ P+	A+ T- P+	−0.068	0.040	−1.706	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A+ T- P+	−0.096	0.043	−2.248	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A+ T- P+	−0.005	0.041	−0.134	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A- T+ P-	−0.189	0.041	−4.604	<0.001	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T+ P-	−0.217	0.044	−4.957	<0.001	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T+ P-	−0.126	0.042	−3.006	0.053	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T+ P+	−0.028	0.043	−0.649	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T+ P+	0.063	0.041	1.521	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A+ T+ P-	0.091	0.044	2.063	≥0.25	Treatment
Richness
Treatment	Alpha.Score	A- T- P+	A- T- P-	0.009	0.034	0.265	≥0.25	Treatment
Treatment	Alpha.Score	A+ T- P-	A- T- P-	−0.049	0.035	−1.430	≥0.25	Treatment
Treatment	Alpha.Score	A+ T- P+	A- T- P-	−0.045	0.034	−1.325	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P-	A- T- P-	−0.016	0.035	−0.467	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A- T- P-	0.016	0.034	0.456	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T- P-	0.029	0.037	0.783	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T- P-	−0.054	0.035	−1.533	≥0.25	Treatment
Treatment	Alpha.Score	A+ T- P-	A- T- P+	−0.058	0.035	−1.691	≥0.25	Treatment
Treatment	Alpha.Score	A+ T- P+	A- T- P+	−0.054	0.034	−1.589	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P-	A- T- P+	−0.025	0.035	−0.724	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A- T- P+	0.007	0.034	0.194	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T- P+	0.020	0.037	0.538	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T- P+	−0.063	0.035	−1.789	≥0.25	Treatment
Treatment	Alpha.Score	A+ T- P+	A+ T- P-	0.004	0.035	0.122	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P-	A+ T- P-	0.033	0.036	0.925	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A+ T- P-	0.065	0.035	1.860	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A+ T- P-	0.078	0.037	2.093	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A+ T- P-	−0.005	0.036	−0.126	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P-	A+ T- P+	0.029	0.035	0.816	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A+ T- P+	0.061	0.035	1.761	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A+ T- P+	0.074	0.037	2.002	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A+ T- P+	−0.009	0.035	−0.248	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A- T+ P-	0.032	0.035	0.903	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T+ P-	0.045	0.038	1.193	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T+ P-	−0.038	0.036	−1.034	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T+ P+	0.013	0.037	0.353	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T+ P+	−0.070	0.036	−1.955	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A+ T+ P-	−0.083	0.038	−2.180	≥0.25	Treatment

60 DPE

Displacement

Plot

Code

# Set the subsection for this analysis
tmp.resSubSection <- "All_60DPE"

# GLM Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]] <- 
  test.alpha.all_metrics$displacement %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models(formula_str = "Displacement ~ Treatment", 
                   alpha_metric_col = "Alpha.Metric", 
                   alpha_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["ANOVA"]] <-
  run_glm_anova(alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]])

# Tukey Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["Tukey"]] <-
  test.alpha.all_metrics$displacement %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm(., "Displacement", "Alpha.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

Tables

GLM Results - 60 DPE
glm(Displacement ~ Treatment); All treatments at 60 DPE
term	estimate	std.error	statistic	p.value	p.adj.sig
Shannon
(Intercept)	0.091	0.037	2.457	0.015	*
A- T- P+	−0.082	0.052	−1.590	0.113	ns
A+ T- P-	−0.034	0.055	−0.617	≥0.25	ns
A+ T- P+	−0.083	0.052	−1.590	0.113	ns
A- T+ P-	−0.041	0.057	−0.719	≥0.25	ns
A- T+ P+	−0.159	0.054	−2.965	0.003	**
A+ T+ P-	−0.191	0.068	−2.815	0.005	**
A+ T+ P+	−0.079	0.058	−1.364	0.174	ns
Simpson
(Intercept)	0.008	0.049	0.162	≥0.25	ns
A- T- P+	−0.070	0.069	−1.011	≥0.25	ns
A+ T- P-	0.023	0.073	0.311	≥0.25	ns
A+ T- P+	−0.070	0.070	−1.002	≥0.25	ns
A- T+ P-	0.030	0.076	0.399	≥0.25	ns
A- T+ P+	−0.249	0.072	−3.483	<0.001	***
A+ T+ P-	−0.332	0.091	−3.665	<0.001	***
A+ T+ P+	−0.106	0.077	−1.387	0.167	ns
Richness
(Intercept)	0.338	0.032	10.566	<0.001	****
A- T- P+	0.025	0.045	0.552	≥0.25	ns
A+ T- P-	−0.064	0.047	−1.349	0.179	ns
A+ T- P+	−0.108	0.046	−2.363	0.019	*
A- T+ P-	−0.089	0.049	−1.804	0.073	ns
A- T+ P+	0.022	0.047	0.462	≥0.25	ns
A+ T+ P-	0.067	0.059	1.141	≥0.25	ns
A+ T+ P+	−0.028	0.050	−0.560	≥0.25	ns

ANOVA of GLM - 60 DPE
ANOVA(GLM(Displacement ~ Treatment), type = 2); All treatments at 60 DPE
term	statistic	df	p.value	sig
Shannon
Treatment	14.694	7.000	0.040	*
Simpson
Treatment	31.743	7.000	<0.001	****
Richness
Treatment	19.090	7.000	0.008	**

Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); All treatments at 60 DPE
term	.y.	group1	group2	estimate	std.error	statistic	adj.p.value	Variable
Shannon
Treatment	Alpha.Score	A- T- P+	A- T- P-	−0.082	0.052	−1.590	≥0.25	Treatment
Treatment	Alpha.Score	A+ T- P-	A- T- P-	−0.034	0.055	−0.617	≥0.25	Treatment
Treatment	Alpha.Score	A+ T- P+	A- T- P-	−0.083	0.052	−1.590	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P-	A- T- P-	−0.041	0.057	−0.719	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A- T- P-	−0.159	0.054	−2.965	0.059	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T- P-	−0.191	0.068	−2.815	0.089	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T- P-	−0.079	0.058	−1.364	≥0.25	Treatment
Treatment	Alpha.Score	A+ T- P-	A- T- P+	0.049	0.054	0.894	≥0.25	Treatment
Treatment	Alpha.Score	A+ T- P+	A- T- P+	−0.001	0.052	−0.022	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P-	A- T- P+	0.041	0.057	0.731	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A- T- P+	−0.077	0.053	−1.442	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T- P+	−0.109	0.068	−1.611	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T- P+	0.004	0.057	0.066	≥0.25	Treatment
Treatment	Alpha.Score	A+ T- P+	A+ T- P-	−0.050	0.055	−0.904	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P-	A+ T- P-	−0.007	0.059	−0.122	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A+ T- P-	−0.126	0.056	−2.234	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A+ T- P-	−0.158	0.070	−2.254	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A+ T- P-	−0.045	0.060	−0.748	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P-	A+ T- P+	0.043	0.057	0.743	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A+ T- P+	−0.076	0.054	−1.402	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A+ T- P+	−0.108	0.068	−1.581	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A+ T- P+	0.005	0.058	0.085	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A- T+ P-	−0.118	0.058	−2.027	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T+ P-	−0.150	0.072	−2.097	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T+ P-	−0.038	0.062	−0.607	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T+ P+	−0.032	0.069	−0.463	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T+ P+	0.081	0.059	1.368	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A+ T+ P-	0.113	0.072	1.562	≥0.25	Treatment
Simpson
Treatment	Alpha.Score	A- T- P+	A- T- P-	−0.070	0.069	−1.011	≥0.25	Treatment
Treatment	Alpha.Score	A+ T- P-	A- T- P-	0.023	0.073	0.311	≥0.25	Treatment
Treatment	Alpha.Score	A+ T- P+	A- T- P-	−0.070	0.070	−1.002	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P-	A- T- P-	0.030	0.076	0.399	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A- T- P-	−0.249	0.072	−3.483	0.012	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T- P-	−0.332	0.091	−3.665	0.006	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T- P-	−0.106	0.077	−1.387	≥0.25	Treatment
Treatment	Alpha.Score	A+ T- P-	A- T- P+	0.092	0.072	1.276	≥0.25	Treatment
Treatment	Alpha.Score	A+ T- P+	A- T- P+	0.000	0.070	−0.005	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P-	A- T- P+	0.100	0.075	1.326	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A- T- P+	−0.180	0.071	−2.525	0.182	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T- P+	−0.262	0.090	−2.906	0.070	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T- P+	−0.037	0.076	−0.480	≥0.25	Treatment
Treatment	Alpha.Score	A+ T- P+	A+ T- P-	−0.093	0.073	−1.265	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P-	A+ T- P-	0.008	0.079	0.096	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A+ T- P-	−0.272	0.075	−3.632	0.007	Treatment
Treatment	Alpha.Score	A+ T+ P-	A+ T- P-	−0.355	0.093	−3.804	0.004	Treatment
Treatment	Alpha.Score	A+ T+ P+	A+ T- P-	−0.129	0.080	−1.617	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P-	A+ T- P+	0.100	0.076	1.315	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A+ T- P+	−0.179	0.072	−2.488	0.198	Treatment
Treatment	Alpha.Score	A+ T+ P-	A+ T- P+	−0.262	0.091	−2.879	0.075	Treatment
Treatment	Alpha.Score	A+ T+ P+	A+ T- P+	−0.036	0.077	−0.470	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A- T+ P-	−0.280	0.078	−3.594	0.007	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T+ P-	−0.362	0.096	−3.790	0.004	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T+ P-	−0.137	0.083	−1.655	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T+ P+	−0.083	0.092	−0.895	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T+ P+	0.143	0.079	1.818	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A+ T+ P-	0.226	0.096	2.343	≥0.25	Treatment
Richness
Treatment	Alpha.Score	A- T- P+	A- T- P-	0.025	0.045	0.552	≥0.25	Treatment
Treatment	Alpha.Score	A+ T- P-	A- T- P-	−0.064	0.047	−1.349	≥0.25	Treatment
Treatment	Alpha.Score	A+ T- P+	A- T- P-	−0.108	0.046	−2.363	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P-	A- T- P-	−0.089	0.049	−1.804	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A- T- P-	0.022	0.047	0.462	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T- P-	0.067	0.059	1.141	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T- P-	−0.028	0.050	−0.560	≥0.25	Treatment
Treatment	Alpha.Score	A+ T- P-	A- T- P+	−0.089	0.047	−1.884	≥0.25	Treatment
Treatment	Alpha.Score	A+ T- P+	A- T- P+	−0.132	0.045	−2.927	0.066	Treatment
Treatment	Alpha.Score	A- T+ P-	A- T- P+	−0.114	0.049	−2.320	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A- T- P+	−0.003	0.046	−0.070	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T- P+	0.043	0.059	0.724	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T- P+	−0.053	0.050	−1.062	≥0.25	Treatment
Treatment	Alpha.Score	A+ T- P+	A+ T- P-	−0.044	0.048	−0.913	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P-	A+ T- P-	−0.025	0.051	−0.488	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A+ T- P-	0.086	0.049	1.755	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A+ T- P-	0.131	0.061	2.164	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A+ T- P-	0.036	0.052	0.694	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P-	A+ T- P+	0.019	0.050	0.373	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A+ T- P+	0.129	0.047	2.752	0.106	Treatment
Treatment	Alpha.Score	A+ T+ P-	A+ T- P+	0.175	0.059	2.954	0.061	Treatment
Treatment	Alpha.Score	A+ T+ P+	A+ T- P+	0.080	0.050	1.585	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A- T+ P-	0.111	0.051	2.184	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T+ P-	0.156	0.062	2.514	0.187	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T+ P-	0.061	0.054	1.137	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T+ P+	0.046	0.060	0.762	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T+ P+	−0.050	0.051	−0.967	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A+ T+ P-	−0.095	0.063	−1.520	≥0.25	Treatment

Significant Results:Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); All treatments at 60 DPE
term	.y.	group1	group2	estimate	std.error	statistic	adj.p.value	Variable
Simpson
Treatment	Alpha.Score	A- T+ P+	A- T- P-	−0.249	0.072	−3.483	0.012	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T- P-	−0.332	0.091	−3.665	0.006	Treatment
Treatment	Alpha.Score	A- T+ P+	A+ T- P-	−0.272	0.075	−3.632	0.007	Treatment
Treatment	Alpha.Score	A+ T+ P-	A+ T- P-	−0.355	0.093	−3.804	0.004	Treatment
Treatment	Alpha.Score	A- T+ P+	A- T+ P-	−0.280	0.078	−3.594	0.007	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T+ P-	−0.362	0.096	−3.790	0.004	Treatment

Beta

All

Displacement

Plot

Code

# Set the subsection for this analysis
tmp.resSubSection <- "All"

# GLM Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]] <- 
  test.beta.all_metrics$volatility %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models_beta(formula_str = "Displacement ~ Treatment", 
                   beta_metric_col = "Beta.Metric", 
                   beta_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["ANOVA"]] <-
  run_glm_anova_beta(beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]])

# Tukey Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["Tukey"]] <-
  test.beta.all_metrics$volatility %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm_beta(., "Displacement", "Beta.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

Tables

GLM Results - Beta
glm(Displacement ~ Treatment); All treatments
term	estimate	std.error	statistic	p.value	p.adj.sig
Bray-Curtis_norm
(Intercept)	−0.061	0.027	−2.266	0.024	*
A- T- P+	0.142	0.038	3.716	<0.001	***
A+ T- P-	0.016	0.039	0.412	≥0.25	ns
A+ T- P+	0.002	0.039	0.040	≥0.25	ns
A- T+ P-	0.081	0.040	2.048	0.041	*
A- T+ P+	−0.163	0.039	−4.216	<0.001	****
A+ T+ P-	−0.020	0.041	−0.473	≥0.25	ns
A+ T+ P+	−0.004	0.040	−0.104	≥0.25	ns
Canberra_norm
(Intercept)	−0.176	0.021	−8.202	<0.001	****
A- T- P+	0.070	0.030	2.300	0.022	*
A+ T- P-	0.040	0.031	1.307	0.192	ns
A+ T- P+	0.083	0.031	2.726	0.007	**
A- T+ P-	0.179	0.031	5.696	<0.001	****
A- T+ P+	0.054	0.031	1.742	0.082	ns
A+ T+ P-	0.155	0.033	4.698	<0.001	****
A+ T+ P+	0.096	0.032	3.044	0.002	**

ANOVA of GLM - Beta
ANOVA(GLM(Displacement ~ Treatment), type = 2); All treatments
term	statistic	df	p.value	sig
Bray-Curtis_norm
Treatment	70.305	7.000	<0.001	****
Canberra_norm
Treatment	46.734	7.000	<0.001	****

Pairwise Tukey's HSD, p.adj: Dunnett - Beta
Tukey(Displacement ~ Treatment); All treatments
term	.y.	group1	group2	estimate	std.error	statistic	adj.p.value	Variable
Bray-Curtis_norm
Treatment	Beta.Score	A- T- P+	A- T- P-	0.142	0.038	3.716	0.005	Treatment
Treatment	Beta.Score	A+ T- P-	A- T- P-	0.016	0.039	0.412	≥0.25	Treatment
Treatment	Beta.Score	A+ T- P+	A- T- P-	0.002	0.039	0.040	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P-	A- T- P-	0.081	0.040	2.048	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P+	A- T- P-	−0.163	0.039	−4.216	<0.001	Treatment
Treatment	Beta.Score	A+ T+ P-	A- T- P-	−0.020	0.041	−0.473	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P+	A- T- P-	−0.004	0.040	−0.104	≥0.25	Treatment
Treatment	Beta.Score	A+ T- P-	A- T- P+	−0.126	0.039	−3.232	0.027	Treatment
Treatment	Beta.Score	A+ T- P+	A- T- P+	−0.141	0.039	−3.652	0.006	Treatment
Treatment	Beta.Score	A- T+ P-	A- T- P+	−0.061	0.040	−1.543	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P+	A- T- P+	−0.306	0.039	−7.884	<0.001	Treatment
Treatment	Beta.Score	A+ T+ P-	A- T- P+	−0.162	0.041	−3.900	0.002	Treatment
Treatment	Beta.Score	A+ T+ P+	A- T- P+	−0.146	0.040	−3.682	0.006	Treatment
Treatment	Beta.Score	A+ T- P+	A+ T- P-	−0.015	0.039	−0.370	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P-	A+ T- P-	0.065	0.040	1.612	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P+	A+ T- P-	−0.180	0.040	−4.543	<0.001	Treatment
Treatment	Beta.Score	A+ T+ P-	A+ T- P-	−0.036	0.042	−0.846	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P+	A+ T- P-	−0.020	0.040	−0.499	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P-	A+ T- P+	0.080	0.040	1.997	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P+	A+ T- P+	−0.165	0.039	−4.230	<0.001	Treatment
Treatment	Beta.Score	A+ T+ P-	A+ T- P+	−0.021	0.042	−0.507	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P+	A+ T- P+	−0.006	0.040	−0.142	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P+	A- T+ P-	−0.245	0.040	−6.102	<0.001	Treatment
Treatment	Beta.Score	A+ T+ P-	A- T+ P-	−0.101	0.043	−2.358	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P+	A- T+ P-	−0.085	0.041	−2.077	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P-	A- T+ P+	0.144	0.042	3.429	0.014	Treatment
Treatment	Beta.Score	A+ T+ P+	A- T+ P+	0.159	0.040	3.961	0.002	Treatment
Treatment	Beta.Score	A+ T+ P+	A+ T+ P-	0.015	0.043	0.361	≥0.25	Treatment
Canberra_norm
Treatment	Beta.Score	A- T- P+	A- T- P-	0.070	0.030	2.300	≥0.25	Treatment
Treatment	Beta.Score	A+ T- P-	A- T- P-	0.040	0.031	1.307	≥0.25	Treatment
Treatment	Beta.Score	A+ T- P+	A- T- P-	0.083	0.031	2.726	0.115	Treatment
Treatment	Beta.Score	A- T+ P-	A- T- P-	0.179	0.031	5.696	<0.001	Treatment
Treatment	Beta.Score	A- T+ P+	A- T- P-	0.054	0.031	1.742	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P-	A- T- P-	0.155	0.033	4.698	<0.001	Treatment
Treatment	Beta.Score	A+ T+ P+	A- T- P-	0.096	0.032	3.044	0.048	Treatment
Treatment	Beta.Score	A+ T- P-	A- T- P+	−0.029	0.031	−0.948	≥0.25	Treatment
Treatment	Beta.Score	A+ T- P+	A- T- P+	0.013	0.031	0.441	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P-	A- T- P+	0.109	0.031	3.473	0.012	Treatment
Treatment	Beta.Score	A- T+ P+	A- T- P+	−0.016	0.031	−0.528	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P-	A- T- P+	0.085	0.033	2.576	0.164	Treatment
Treatment	Beta.Score	A+ T+ P+	A- T- P+	0.026	0.032	0.830	≥0.25	Treatment
Treatment	Beta.Score	A+ T- P+	A+ T- P-	0.043	0.031	1.374	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P-	A+ T- P-	0.138	0.032	4.327	<0.001	Treatment
Treatment	Beta.Score	A- T+ P+	A+ T- P-	0.013	0.031	0.418	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P-	A+ T- P-	0.114	0.033	3.410	0.015	Treatment
Treatment	Beta.Score	A+ T+ P+	A+ T- P-	0.055	0.032	1.729	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P-	A+ T- P+	0.096	0.032	3.027	0.050	Treatment
Treatment	Beta.Score	A- T+ P+	A+ T- P+	−0.030	0.031	−0.960	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P-	A+ T- P+	0.071	0.033	2.156	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P+	A+ T- P+	0.013	0.032	0.400	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P+	A- T+ P-	−0.125	0.032	−3.943	0.002	Treatment
Treatment	Beta.Score	A+ T+ P-	A- T+ P-	−0.024	0.034	−0.718	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P+	A- T+ P-	−0.083	0.033	−2.549	0.175	Treatment
Treatment	Beta.Score	A+ T+ P-	A- T+ P+	0.101	0.033	3.036	0.049	Treatment
Treatment	Beta.Score	A+ T+ P+	A- T+ P+	0.042	0.032	1.329	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P+	A+ T+ P-	−0.059	0.034	−1.725	≥0.25	Treatment

60 DPE

Displacement

Plot

Code

# Set the subsection for this analysis
tmp.resSubSection <- "All_60DPE"

# GLM Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]] <- 
  test.beta.all_metrics$displacement %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models_beta(formula_str = "Displacement ~ Treatment", 
                   beta_metric_col = "Beta.Metric", 
                   beta_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["ANOVA"]] <-
  run_glm_anova_beta(beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]])

# Tukey Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["Tukey"]] <-
  test.beta.all_metrics$displacement %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm_beta(., "Displacement", "Beta.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 

## Warning in RET$pfunction("adjusted", ...): Completion with error > abseps

Tables

GLM Results - 60 DPE
glm(Displacement ~ Treatment); All treatments at 60 DPE
term	estimate	std.error	statistic	p.value	p.adj.sig
Bray-Curtis_norm
(Intercept)	−0.060	0.041	−1.474	0.142	ns
A- T- P+	0.284	0.057	4.945	<0.001	****
A+ T- P-	0.007	0.061	0.115	≥0.25	ns
A+ T- P+	−0.021	0.058	−0.360	≥0.25	ns
A- T+ P-	0.043	0.063	0.675	≥0.25	ns
A- T+ P+	−0.206	0.060	−3.458	<0.001	***
A+ T+ P-	−0.119	0.075	−1.584	0.115	ns
A+ T+ P+	0.063	0.064	0.993	≥0.25	ns
Canberra_norm
(Intercept)	−0.124	0.034	−3.647	<0.001	***
A- T- P+	0.134	0.048	2.796	0.006	**
A+ T- P-	−0.089	0.051	−1.750	0.081	ns
A+ T- P+	0.037	0.049	0.765	≥0.25	ns
A- T+ P-	0.099	0.053	1.871	0.063	ns
A- T+ P+	0.065	0.050	1.306	0.193	ns
A+ T+ P-	0.123	0.063	1.962	0.051	ns
A+ T+ P+	0.071	0.053	1.329	0.185	ns

ANOVA of GLM - 60 DPE
ANOVA(GLM(Displacement ~ Treatment), type = 2); All treatments at 60 DPE
term	statistic	df	p.value	sig
Bray-Curtis_norm
Treatment	77.159	7.000	<0.001	****
Canberra_norm
Treatment	26.083	7.000	<0.001	***

Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); All treatments at 60 DPE
term	.y.	group1	group2	estimate	std.error	statistic	adj.p.value	Variable
Bray-Curtis_norm
Treatment	Beta.Score	A- T- P+	A- T- P-	0.284	0.057	4.945	<0.001	Treatment
Treatment	Beta.Score	A+ T- P-	A- T- P-	0.007	0.061	0.115	≥0.25	Treatment
Treatment	Beta.Score	A+ T- P+	A- T- P-	−0.021	0.058	−0.360	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P-	A- T- P-	0.043	0.063	0.675	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P+	A- T- P-	−0.206	0.060	−3.458	0.013	Treatment
Treatment	Beta.Score	A+ T+ P-	A- T- P-	−0.119	0.075	−1.584	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P+	A- T- P-	0.063	0.064	0.993	≥0.25	Treatment
Treatment	Beta.Score	A+ T- P-	A- T- P+	−0.277	0.060	−4.595	<0.001	Treatment
Treatment	Beta.Score	A+ T- P+	A- T- P+	−0.305	0.058	−5.272	<0.001	Treatment
Treatment	Beta.Score	A- T+ P-	A- T- P+	−0.241	0.063	−3.845	0.003	Treatment
Treatment	Beta.Score	A- T+ P+	A- T- P+	−0.490	0.059	−8.275	<0.001	Treatment
Treatment	Beta.Score	A+ T+ P-	A- T- P+	−0.403	0.075	−5.372	<0.001	Treatment
Treatment	Beta.Score	A+ T+ P+	A- T- P+	−0.220	0.063	−3.473	0.012	Treatment
Treatment	Beta.Score	A+ T- P+	A+ T- P-	−0.028	0.061	−0.458	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P-	A+ T- P-	0.036	0.066	0.542	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P+	A+ T- P-	−0.213	0.062	−3.418	0.014	Treatment
Treatment	Beta.Score	A+ T+ P-	A+ T- P-	−0.126	0.078	−1.629	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P+	A+ T- P-	0.056	0.066	0.849	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P-	A+ T- P+	0.064	0.063	1.001	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P+	A+ T- P+	−0.185	0.060	−3.086	0.042	Treatment
Treatment	Beta.Score	A+ T+ P-	A+ T- P+	−0.098	0.076	−1.300	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P+	A+ T- P+	0.084	0.064	1.313	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P+	A- T+ P-	−0.249	0.065	−3.840	0.003	Treatment
Treatment	Beta.Score	A+ T+ P-	A- T+ P-	−0.162	0.079	−2.037	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P+	A- T+ P-	0.021	0.069	0.303	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P-	A- T+ P+	0.087	0.077	1.130	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P+	A- T+ P+	0.269	0.065	4.116	<0.001	Treatment
Treatment	Beta.Score	A+ T+ P+	A+ T+ P-	0.183	0.080	2.282	≥0.25	Treatment
Canberra_norm
Treatment	Beta.Score	A- T- P+	A- T- P-	0.134	0.048	2.796	0.094	Treatment
Treatment	Beta.Score	A+ T- P-	A- T- P-	−0.089	0.051	−1.750	≥0.25	Treatment
Treatment	Beta.Score	A+ T- P+	A- T- P-	0.037	0.049	0.765	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P-	A- T- P-	0.099	0.053	1.871	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P+	A- T- P-	0.065	0.050	1.306	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P-	A- T- P-	0.123	0.063	1.962	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P+	A- T- P-	0.071	0.053	1.329	≥0.25	Treatment
Treatment	Beta.Score	A+ T- P-	A- T- P+	−0.223	0.050	−4.425	<0.001	Treatment
Treatment	Beta.Score	A+ T- P+	A- T- P+	−0.097	0.048	−2.006	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P-	A- T- P+	−0.035	0.052	−0.676	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P+	A- T- P+	−0.069	0.049	−1.397	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P-	A- T- P+	−0.011	0.063	−0.168	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P+	A- T- P+	−0.063	0.053	−1.192	≥0.25	Treatment
Treatment	Beta.Score	A+ T- P+	A+ T- P-	0.126	0.051	2.468	0.206	Treatment
Treatment	Beta.Score	A- T+ P-	A+ T- P-	0.187	0.055	3.412	0.014	Treatment
Treatment	Beta.Score	A- T+ P+	A+ T- P-	0.153	0.052	2.951	0.062	Treatment
Treatment	Beta.Score	A+ T+ P-	A+ T- P-	0.212	0.065	3.275	0.023	Treatment
Treatment	Beta.Score	A+ T+ P+	A+ T- P-	0.159	0.055	2.875	0.076	Treatment
Treatment	Beta.Score	A- T+ P-	A+ T- P+	0.061	0.053	1.159	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P+	A+ T- P+	0.028	0.050	0.555	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P-	A+ T- P+	0.086	0.063	1.366	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P+	A+ T- P+	0.034	0.054	0.628	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P+	A- T+ P-	−0.034	0.054	−0.623	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P-	A- T+ P-	0.025	0.066	0.375	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P+	A- T+ P-	−0.028	0.057	−0.485	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P-	A- T+ P+	0.058	0.064	0.913	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P+	A- T+ P+	0.006	0.055	0.107	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P+	A+ T+ P-	−0.053	0.067	−0.788	≥0.25	Treatment

Significant Results:Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); All treatments at 60 DPE
term	.y.	group1	group2	estimate	std.error	statistic	adj.p.value	Variable
Bray-Curtis_norm
Treatment	Beta.Score	A- T- P+	A- T- P-	0.284	0.057	4.945	<0.001	Treatment
Treatment	Beta.Score	A- T+ P+	A- T- P-	−0.206	0.060	−3.458	0.013	Treatment
Treatment	Beta.Score	A+ T- P-	A- T- P+	−0.277	0.060	−4.595	<0.001	Treatment
Treatment	Beta.Score	A+ T- P+	A- T- P+	−0.305	0.058	−5.272	<0.001	Treatment
Treatment	Beta.Score	A- T+ P-	A- T- P+	−0.241	0.063	−3.845	0.003	Treatment
Treatment	Beta.Score	A- T+ P+	A- T- P+	−0.490	0.059	−8.275	<0.001	Treatment
Treatment	Beta.Score	A+ T+ P-	A- T- P+	−0.403	0.075	−5.372	<0.001	Treatment
Treatment	Beta.Score	A+ T+ P+	A- T- P+	−0.220	0.063	−3.473	0.012	Treatment
Treatment	Beta.Score	A- T+ P+	A+ T- P-	−0.213	0.062	−3.418	0.014	Treatment
Treatment	Beta.Score	A- T+ P+	A+ T- P+	−0.185	0.060	−3.086	0.042	Treatment
Treatment	Beta.Score	A- T+ P+	A- T+ P-	−0.249	0.065	−3.840	0.003	Treatment
Treatment	Beta.Score	A+ T+ P+	A- T+ P+	0.269	0.065	4.116	<0.001	Treatment
Canberra_norm
Treatment	Beta.Score	A+ T- P-	A- T- P+	−0.223	0.050	−4.425	<0.001	Treatment
Treatment	Beta.Score	A- T+ P-	A+ T- P-	0.187	0.055	3.412	0.014	Treatment
Treatment	Beta.Score	A+ T+ P-	A+ T- P-	0.212	0.065	3.275	0.023	Treatment

PCoA

Plot

Code

# Set the subsection for this analysis
tmp.resSubSection <- "All_60DPE"



#### Dispersion --------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]] <-
  run_BetaDispersion(dist.matrix = beta.dist.mat[[tmp.resSubSection]], 
                     var = c("treatment_group"))

##### ANOVA ------------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA"]] <-
  get_HoD_anova(betaDisper = beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]],
                var = c("treatment_group"))

### Table

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA.Table"]] <-
  beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA"]] %>%
  
    # Create the table
    dplyr::group_by(Beta.Metric) %>%
    set_GT(var = "p.value", group.by = "Beta.Metric")  %>%
  
    # Title/caption
    gt::tab_header(
      title = "ANOVA: Homogeneity of Dispersion",
      subtitle = "ANOVA(Beta Disperson ~ Treatment); All treatments at 60 DPE"
    ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = c("p.value", "sig"),
      rows = p.value < 0.05
    )
  )

##### Tukey ------------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey"]] <-
  get_HoD_tukey(betaDisper = beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]],
                var = c("treatment_group")) 

### Table

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey.Table"]] <-
  beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey"]] %>%
    
    dplyr::mutate(
        group1 = case_when(
        group1 == "Control" ~ "A- T- P-",
        group1 == "Parasite" ~ "A- T- P+",
        group1 == "Antibiotics" ~ "A+ T- P-",
        group1 == "Antibiotics_Parasite" ~ "A+ T- P+",
        group1 == "Temperature" ~ "A- T+ P-",
        group1 == "Temperature_Parasite" ~ "A- T+ P+",
        group1 == "Antibiotics_Temperature" ~ "A+ T+ P-",
        group1 == "Antibiotics_Temperature_Parasite" ~ "A+ T+ P+",
    ),
        group2 = case_when(
        group2 == "Control" ~ "A- T- P-",
        group2 == "Parasite" ~ "A- T- P+",
        group2 == "Antibiotics" ~ "A+ T- P-",
        group2 == "Antibiotics_Parasite" ~ "A+ T- P+",
        group2 == "Temperature" ~ "A- T+ P-",
        group2 == "Temperature_Parasite" ~ "A- T+ P+",
        group2 == "Antibiotics_Temperature" ~ "A+ T+ P-",
        group2 == "Antibiotics_Temperature_Parasite" ~ "A+ T+ P+",
    )) %>%
  
  # Create the table
  dplyr::group_by(Beta.Metric) %>%
  set_GT(var = "adj.p.value", group.by = "Beta.Metric")  %>%
  
  # Title/caption
  gt::tab_header(
    title = "Tukey: Homogeneity of Dispersion",
    subtitle = "Tukey(Beta Disperson ~ Treatment); All treatments at 60 DPE"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = c("adj.p.value", "sig"),
      rows = adj.p.value < 0.05
    )
  )

Tables

ANOVA: Homogeneity of Dispersion
ANOVA(Beta Disperson ~ Treatment); All treatments at 60 DPE
term	df	sumsq	meansq	statistic	p.value	sig
bray
treatment_group	7.000	0.645	0.092	3.523	0.001	**
Residual	228.000	5.965	0.026	NA	NA	NA
canberra
treatment_group	7.000	0.132	0.019	5.021	<0.001	****
Residual	228.000	0.857	0.004	NA	NA	NA

Tukey: Homogeneity of Dispersion
Tukey(Beta Disperson ~ Treatment); All treatments at 60 DPE
.y.	term	group1	group2	estimate	conf.low	conf.high	adj.p.value	sig
bray
Distance	treatment_group	A- T- P+	A- T- P-	−0.025	−0.141	0.091	≥0.25	ns
Distance	treatment_group	A+ T- P-	A- T- P-	−0.036	−0.158	0.086	≥0.25	ns
Distance	treatment_group	A+ T- P+	A- T- P-	−0.065	−0.183	0.052	≥0.25	ns
Distance	treatment_group	A- T+ P-	A- T- P-	−0.064	−0.191	0.064	≥0.25	ns
Distance	treatment_group	A- T+ P+	A- T- P-	−0.106	−0.227	0.014	0.126	ns
Distance	treatment_group	A+ T+ P-	A- T- P-	−0.162	−0.314	−0.010	0.028	*
Distance	treatment_group	A+ T+ P+	A- T- P-	−0.157	−0.286	−0.028	0.006	**
Distance	treatment_group	A+ T- P-	A- T- P+	−0.011	−0.132	0.111	≥0.25	ns
Distance	treatment_group	A+ T- P+	A- T- P+	−0.040	−0.157	0.076	≥0.25	ns
Distance	treatment_group	A- T+ P-	A- T- P+	−0.039	−0.165	0.088	≥0.25	ns
Distance	treatment_group	A- T+ P+	A- T- P+	−0.081	−0.201	0.038	≥0.25	ns
Distance	treatment_group	A+ T+ P-	A- T- P+	−0.137	−0.289	0.014	0.108	ns
Distance	treatment_group	A+ T+ P+	A- T- P+	−0.132	−0.260	−0.004	0.038	*
Distance	treatment_group	A+ T- P+	A+ T- P-	−0.029	−0.152	0.094	≥0.25	ns
Distance	treatment_group	A- T+ P-	A+ T- P-	−0.028	−0.160	0.105	≥0.25	ns
Distance	treatment_group	A- T+ P+	A+ T- P-	−0.070	−0.196	0.055	≥0.25	ns
Distance	treatment_group	A+ T+ P-	A+ T- P-	−0.126	−0.283	0.030	0.215	ns
Distance	treatment_group	A+ T+ P+	A+ T- P-	−0.121	−0.255	0.013	0.108	ns
Distance	treatment_group	A- T+ P-	A+ T- P+	0.001	−0.127	0.130	≥0.25	ns
Distance	treatment_group	A- T+ P+	A+ T- P+	−0.041	−0.162	0.080	≥0.25	ns
Distance	treatment_group	A+ T+ P-	A+ T- P+	−0.097	−0.250	0.056	≥0.25	ns
Distance	treatment_group	A+ T+ P+	A+ T- P+	−0.092	−0.222	0.038	≥0.25	ns
Distance	treatment_group	A- T+ P+	A- T+ P-	−0.042	−0.173	0.088	≥0.25	ns
Distance	treatment_group	A+ T+ P-	A- T+ P-	−0.098	−0.259	0.062	≥0.25	ns
Distance	treatment_group	A+ T+ P+	A- T+ P-	−0.093	−0.232	0.045	≥0.25	ns
Distance	treatment_group	A+ T+ P-	A- T+ P+	−0.056	−0.211	0.099	≥0.25	ns
Distance	treatment_group	A+ T+ P+	A- T+ P+	−0.051	−0.183	0.081	≥0.25	ns
Distance	treatment_group	A+ T+ P+	A+ T+ P-	0.005	−0.157	0.166	≥0.25	ns
canberra
Distance	treatment_group	A- T- P+	A- T- P-	0.008	−0.036	0.052	≥0.25	ns
Distance	treatment_group	A+ T- P-	A- T- P-	−0.056	−0.103	−0.010	0.006	**
Distance	treatment_group	A+ T- P+	A- T- P-	0.007	−0.037	0.052	≥0.25	ns
Distance	treatment_group	A- T+ P-	A- T- P-	0.010	−0.039	0.058	≥0.25	ns
Distance	treatment_group	A- T+ P+	A- T- P-	−0.007	−0.052	0.039	≥0.25	ns
Distance	treatment_group	A+ T+ P-	A- T- P-	−0.023	−0.080	0.035	≥0.25	ns
Distance	treatment_group	A+ T+ P+	A- T- P-	−0.044	−0.093	0.005	0.116	ns
Distance	treatment_group	A+ T- P-	A- T- P+	−0.065	−0.111	−0.019	<0.001	***
Distance	treatment_group	A+ T- P+	A- T- P+	−0.001	−0.045	0.043	≥0.25	ns
Distance	treatment_group	A- T+ P-	A- T- P+	0.001	−0.047	0.049	≥0.25	ns
Distance	treatment_group	A- T+ P+	A- T- P+	−0.015	−0.060	0.030	≥0.25	ns
Distance	treatment_group	A+ T+ P-	A- T- P+	−0.031	−0.088	0.026	≥0.25	ns
Distance	treatment_group	A+ T+ P+	A- T- P+	−0.052	−0.101	−0.004	0.026	*
Distance	treatment_group	A+ T- P+	A+ T- P-	0.064	0.017	0.110	0.001	**
Distance	treatment_group	A- T+ P-	A+ T- P-	0.066	0.016	0.116	0.002	**
Distance	treatment_group	A- T+ P+	A+ T- P-	0.050	0.002	0.097	0.034	*
Distance	treatment_group	A+ T+ P-	A+ T- P-	0.034	−0.026	0.093	≥0.25	ns
Distance	treatment_group	A+ T+ P+	A+ T- P-	0.013	−0.038	0.063	≥0.25	ns
Distance	treatment_group	A- T+ P-	A+ T- P+	0.002	−0.046	0.051	≥0.25	ns
Distance	treatment_group	A- T+ P+	A+ T- P+	−0.014	−0.060	0.032	≥0.25	ns
Distance	treatment_group	A+ T+ P-	A+ T- P+	−0.030	−0.088	0.028	≥0.25	ns
Distance	treatment_group	A+ T+ P+	A+ T- P+	−0.051	−0.100	−0.002	0.035	*
Distance	treatment_group	A- T+ P+	A- T+ P-	−0.016	−0.066	0.033	≥0.25	ns
Distance	treatment_group	A+ T+ P-	A- T+ P-	−0.032	−0.093	0.029	≥0.25	ns
Distance	treatment_group	A+ T+ P+	A- T+ P-	−0.053	−0.106	−0.001	0.044	*
Distance	treatment_group	A+ T+ P-	A- T+ P+	−0.016	−0.075	0.043	≥0.25	ns
Distance	treatment_group	A+ T+ P+	A- T+ P+	−0.037	−0.087	0.013	≥0.25	ns
Distance	treatment_group	A+ T+ P+	A+ T+ P-	−0.021	−0.082	0.040	≥0.25	ns

CAP

Plot

Code

# Set the subsection for this analysis
tmp.resSubSection <- "All_60DPE"


#### Capscale ----------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.mod"]] <-
  run_capscale(ps.cleaned_alpha %>% # is the phyloseq object we are looping over
          # microViz::ps_filter(Treatment  %in% c(
          # "A- T- P-",
          # "A- T- P+",
          # "A+ T- P-",
          # "A+ T- P+",
          # "A- T+ P-",
          # "A- T+ P+",
          # "A+ T+ P-",
          # "A+ T+ P+"
          # )) %>%
          microViz::ps_filter(Time == 60), 
               dist.matrix = beta.dist.mat[[tmp.resSubSection]], 
               formula_str = "dist ~ Treatment")

##### ADONIS ------------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS"]] <-
  run_cap_adonis(ps.cleaned_alpha %>% # is the phyloseq object we are looping over
          # microViz::ps_filter(Treatment  %in% c(
          # "A- T- P-",
          # "A- T- P+",
          # "A+ T- P-",
          # "A+ T- P+",
          # "A- T+ P-",
          # "A- T+ P+",
          # "A+ T+ P-",
          # "A+ T+ P+"
          # )) %>%
          microViz::ps_filter(Time == 60),
                 dist.matrix = beta.dist.mat[[tmp.resSubSection]], 
                 formula_str = "dist ~ Treatment",
                 by.method = "terms") 

## Warning in tidy.anova(.): The following column names in ANOVA output were not
## recognized or transformed: SumOfSqs, R2

## Warning in tidy.anova(.): The following column names in ANOVA output were not
## recognized or transformed: SumOfSqs, R2

### Table

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS.Table"]] <-
  beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS"]] %>%
    set_GT(var = "p.value", group.by = "Beta.Metric") %>%
  
  # Title/caption
  gt::tab_header(
    title = "ADONIS2",
    subtitle = "adonis2(Beta Distance ~ Treatment); All treatments at 60 DPE"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = c("p.value", "sig"),
      rows = p.value < 0.05
    )
  )

Tables

ADONIS2
adonis2(Beta Distance ~ Treatment); All treatments at 60 DPE
term	df	SumOfSqs	R2	statistic	p.value	sig
bray
Treatment	7.000	3.508	0.105	3.813	0.001	***
Residual	228.000	29.963	0.895	NA	NA	NA
Total	235.000	33.471	1.000	NA	NA	NA
canberra
Treatment	7.000	8.648	0.103	3.733	0.001	***
Residual	228.000	75.464	0.897	NA	NA	NA
Total	235.000	84.112	1.000	NA	NA	NA

01) What is the host-microbiome response to parasite exposure?

Alpha

All

Displacement

Plot

Code

# Set the subsection for this analysis
tmp.resSubSection <- "ParasiteExp"

# GLM Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]] <- 
  test.alpha.all_metrics$volatility %>%
    # Filter 
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          "A- T- P+"#,
          # "A+ T- P-",
          # "A+ T- P+",
          # "A- T+ P-",
          # "A- T+ P+",
          # "A+ T+ P-",
          # "A+ T+ P+"
          )) %>%
    
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models(formula_str = "Displacement ~ Treatment", 
                   alpha_metric_col = "Alpha.Metric", 
                   alpha_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["ANOVA"]] <-
  run_glm_anova(alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]])

# Tukey Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["Tukey"]] <-
  test.alpha.all_metrics$volatility %>%
    # Filter 
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          "A- T- P+"#,
          # "A+ T- P-",
          # "A+ T- P+",
          # "A- T+ P-",
          # "A- T+ P+",
          # "A+ T+ P-",
          # "A+ T+ P+"
          )) %>%
    
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm(., "Displacement", "Alpha.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 

Tables

GLM Results
glm(Displacement ~ Treatment); Parasite exposure only
term	estimate	std.error	statistic	p.value	p.adj.sig
Shannon
(Intercept)	0.078	0.021	3.755	<0.001	***
A- T- P+	−0.045	0.029	−1.529	0.128	ns
Simpson
(Intercept)	0.042	0.029	1.480	0.141	ns
A- T- P+	−0.049	0.040	−1.199	0.232	ns
Richness
(Intercept)	0.221	0.025	8.745	<0.001	****
A- T- P+	0.009	0.036	0.251	≥0.25	ns

ANOVA of GLM
ANOVA(GLM(Displacement ~ Treatment), type = 2); Parasite exposure only
term	statistic	df	p.value	sig
Shannon
Treatment	2.338	1.000	0.126	ns
Simpson
Treatment	1.439	1.000	0.230	ns
Richness
Treatment	0.063	1.000	≥0.25	ns

Pairwise Tukey's HSD, p.adj: Dunnett
Tukey(Displacement ~ Treatment); Parasite exposure only
term	.y.	group1	group2	estimate	std.error	statistic	adj.p.value	Variable
Shannon
Treatment	Alpha.Score	A- T- P+	A- T- P-	−0.045	0.029	−1.529	0.126	Treatment
Simpson
Treatment	Alpha.Score	A- T- P+	A- T- P-	−0.049	0.040	−1.199	0.230	Treatment
Richness
Treatment	Alpha.Score	A- T- P+	A- T- P-	0.009	0.036	0.251	≥0.25	Treatment

60 DPE

Displacement

Plot

Code

# Set the subsection for this analysis
tmp.resSubSection <- "ParasiteExp_60DPE"

# GLM Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]] <- 
  test.alpha.all_metrics$displacement %>%
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          "A- T- P+"#,
          # "A+ T- P-",
          # "A+ T- P+",
          # "A- T+ P-",
          # "A- T+ P+",
          # "A+ T+ P-",
          # "A+ T+ P+"
          )) %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models(formula_str = "Displacement ~ Treatment", 
                   alpha_metric_col = "Alpha.Metric", 
                   alpha_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["ANOVA"]] <-
  run_glm_anova(alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]])

# Tukey Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["Tukey"]] <-
  test.alpha.all_metrics$displacement %>%
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          "A- T- P+"#,
          # "A+ T- P-",
          # "A+ T- P+",
          # "A- T+ P-",
          # "A- T+ P+",
          # "A+ T+ P-",
          # "A+ T+ P+"
          )) %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm(., "Displacement", "Alpha.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 

Tables

GLM Results - 60 DPE
glm(Displacement ~ Treatment); All treatments at 60 DPE
term	estimate	std.error	statistic	p.value	p.adj.sig
Shannon
(Intercept)	0.091	0.041	2.207	0.031	*
A- T- P+	−0.082	0.058	−1.428	0.158	ns
Simpson
(Intercept)	0.008	0.053	0.150	≥0.25	ns
A- T- P+	−0.070	0.075	−0.933	≥0.25	ns
Richness
(Intercept)	0.338	0.034	9.916	<0.001	****
A- T- P+	0.025	0.048	0.518	≥0.25	ns

ANOVA of GLM - 60 DPE
ANOVA(GLM(Displacement ~ Treatment), type = 2); Parasite exposure only at 60 DPE
term	statistic	df	p.value	sig
Shannon
Treatment	2.039	1.000	0.153	ns
Simpson
Treatment	0.871	1.000	≥0.25	ns
Richness
Treatment	0.269	1.000	≥0.25	ns

Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); Parasite exposure only at 60 DPE
term	.y.	group1	group2	estimate	std.error	statistic	adj.p.value	Variable
Shannon
Treatment	Alpha.Score	A- T- P+	A- T- P-	−0.082	0.058	−1.428	0.153	Treatment
Simpson
Treatment	Alpha.Score	A- T- P+	A- T- P-	−0.070	0.075	−0.933	≥0.25	Treatment
Richness
Treatment	Alpha.Score	A- T- P+	A- T- P-	0.025	0.048	0.518	≥0.25	Treatment

Significant Results:Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); All treatments at 60 DPE
term	.y.	group1	group2	estimate	std.error	statistic	adj.p.value	Variable

Beta

All

Displacement

Plot

Code

# Set the subsection for this analysis
tmp.resSubSection <- "ParasiteExp"

# GLM Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]] <- 
  test.beta.all_metrics$volatility %>%
        # Filter 
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          "A- T- P+"#,
          # "A+ T- P-",
          # "A+ T- P+",
          # "A- T+ P-",
          # "A- T+ P+",
          # "A+ T+ P-",
          # "A+ T+ P+"
          )) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models_beta(formula_str = "Displacement ~ Treatment", 
                   beta_metric_col = "Beta.Metric", 
                   beta_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["ANOVA"]] <-
  run_glm_anova_beta(beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]])

# Tukey Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["Tukey"]] <-
  test.beta.all_metrics$volatility %>%
        # Filter 
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          "A- T- P+"#,
          # "A+ T- P-",
          # "A+ T- P+",
          # "A- T+ P-",
          # "A- T+ P+",
          # "A+ T+ P-",
          # "A+ T+ P+"
          )) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm_beta(., "Displacement", "Beta.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 

Tables

GLM Results - Beta
glm(Displacement ~ Treatment); Parasite exposure only
term	estimate	std.error	statistic	p.value	p.adj.sig
Bray-Curtis_norm
(Intercept)	−0.061	0.026	−2.369	0.019	*
A- T- P+	0.142	0.037	3.884	<0.001	***
Canberra_norm
(Intercept)	−0.176	0.023	−7.710	<0.001	****
A- T- P+	0.070	0.032	2.162	0.032	*

ANOVA of GLM - Beta
ANOVA(GLM(Displacement ~ Treatment), type = 2); Parasite exposure only
term	statistic	df	p.value	sig
Bray-Curtis_norm
Treatment	15.083	1.000	<0.001	***
Canberra_norm
Treatment	4.673	1.000	0.031	*

Pairwise Tukey's HSD, p.adj: Dunnett - Beta
Tukey(Displacement ~ Treatment); Parasite exposure only
term	.y.	group1	group2	estimate	std.error	statistic	adj.p.value	Variable
Bray-Curtis_norm
Treatment	Beta.Score	A- T- P+	A- T- P-	0.142	0.037	3.884	<0.001	Treatment
Canberra_norm
Treatment	Beta.Score	A- T- P+	A- T- P-	0.070	0.032	2.162	0.031	Treatment

60 DPE

Displacement

Plot

Code

# Set the subsection for this analysis
tmp.resSubSection <- "ParasiteExp_60DPE"

# GLM Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]] <- 
  test.beta.all_metrics$displacement %>%
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          "A- T- P+"#,
          # "A+ T- P-",
          # "A+ T- P+",
          # "A- T+ P-",
          # "A- T+ P+",
          # "A+ T+ P-",
          # "A+ T+ P+"
          )) %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models_beta(formula_str = "Displacement ~ Treatment", 
                   beta_metric_col = "Beta.Metric", 
                   beta_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["ANOVA"]] <-
  run_glm_anova_beta(beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]])

# Tukey Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["Tukey"]] <-
  test.beta.all_metrics$displacement %>%
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          "A- T- P+"#,
          # "A+ T- P-",
          # "A+ T- P+",
          # "A- T+ P-",
          # "A- T+ P+",
          # "A+ T+ P-",
          # "A+ T+ P+"
          )) %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm_beta(., "Displacement", "Beta.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 

Tables

GLM Results - 60 DPE
glm(Displacement ~ Treatment); Parasite exposure only at 60 DPE
term	estimate	std.error	statistic	p.value	p.adj.sig
Bray-Curtis_norm
(Intercept)	−0.060	0.036	−1.654	0.103	ns
A- T- P+	0.284	0.051	5.548	<0.001	****
Canberra_norm
(Intercept)	−0.124	0.028	−4.457	<0.001	****
A- T- P+	0.134	0.039	3.417	0.001	**

ANOVA of GLM - 60 DPE
ANOVA(GLM(Displacement ~ Treatment), type = 2); Parasite exposure only at 60 DPE
term	statistic	df	p.value	sig
Bray-Curtis_norm
Treatment	30.779	1.000	<0.001	****
Canberra_norm
Treatment	11.677	1.000	<0.001	***

Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); Parasite exposure only at 60 DPE
term	.y.	group1	group2	estimate	std.error	statistic	adj.p.value	Variable
Bray-Curtis_norm
Treatment	Beta.Score	A- T- P+	A- T- P-	0.284	0.051	5.548	<0.001	Treatment
Canberra_norm
Treatment	Beta.Score	A- T- P+	A- T- P-	0.134	0.039	3.417	<0.001	Treatment

Significant Results:Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); All treatments at 60 DPE
term	.y.	group1	group2	estimate	std.error	statistic	adj.p.value	Variable
Bray-Curtis_norm
Treatment	Beta.Score	A- T- P+	A- T- P-	0.284	0.051	5.548	<0.001	Treatment
Canberra_norm
Treatment	Beta.Score	A- T- P+	A- T- P-	0.134	0.039	3.417	<0.001	Treatment

PCoA

Plot

Code

# Set the subsection for this analysis
tmp.resSubSection <- "ParasiteExp_60DPE"



#### Dispersion --------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]] <-
  run_BetaDispersion(dist.matrix = beta.dist.mat[[tmp.resSubSection]], 
                     var = c("treatment_group"))

##### ANOVA ------------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA"]] <-
  get_HoD_anova(betaDisper = beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]],
                var = c("treatment_group"))

### Table

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA.Table"]] <-
  beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA"]] %>%
  
    # Create the table
    dplyr::group_by(Beta.Metric) %>%
    set_GT(var = "p.value", group.by = "Beta.Metric")  %>%
  
    # Title/caption
    gt::tab_header(
      title = "ANOVA: Homogeneity of Dispersion",
      subtitle = "ANOVA(Beta Disperson ~ Treatment); Parasite exposure only at 60 DPE"
    ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = c("p.value", "sig"),
      rows = p.value < 0.05
    )
  )

##### Tukey ------------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey"]] <-
  get_HoD_tukey(betaDisper = beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]],
                var = c("treatment_group")) 

### Table

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey.Table"]] <-
  beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey"]] %>%
    
    dplyr::mutate(
        group1 = case_when(
        group1 == "Control" ~ "A- T- P-",
        group1 == "Parasite" ~ "A- T- P+",
        group1 == "Antibiotics" ~ "A+ T- P-",
        group1 == "Antibiotics_Parasite" ~ "A+ T- P+",
        group1 == "Temperature" ~ "A- T+ P-",
        group1 == "Temperature_Parasite" ~ "A- T+ P+",
        group1 == "Antibiotics_Temperature" ~ "A+ T+ P-",
        group1 == "Antibiotics_Temperature_Parasite" ~ "A+ T+ P+",
    ),
        group2 = case_when(
        group2 == "Control" ~ "A- T- P-",
        group2 == "Parasite" ~ "A- T- P+",
        group2 == "Antibiotics" ~ "A+ T- P-",
        group2 == "Antibiotics_Parasite" ~ "A+ T- P+",
        group2 == "Temperature" ~ "A- T+ P-",
        group2 == "Temperature_Parasite" ~ "A- T+ P+",
        group2 == "Antibiotics_Temperature" ~ "A+ T+ P-",
        group2 == "Antibiotics_Temperature_Parasite" ~ "A+ T+ P+",
    )) %>%
  
  # Create the table
  dplyr::group_by(Beta.Metric) %>%
  set_GT(var = "adj.p.value", group.by = "Beta.Metric")  %>%
  
  # Title/caption
  gt::tab_header(
    title = "Tukey: Homogeneity of Dispersion",
    subtitle = "Tukey(Beta Disperson ~ Treatment); Parasite exposure only 60 DPE"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = c("adj.p.value", "sig"),
      rows = adj.p.value < 0.05
    )
  )

Tables

ANOVA: Homogeneity of Dispersion
ANOVA(Beta Disperson ~ Treatment); Parasite exposure only at 60 DPE
term	df	sumsq	meansq	statistic	p.value	sig
bray
treatment_group	1.000	0.011	0.011	0.388	≥0.25	ns
Residual	71.000	2.088	0.029	NA	NA	NA
canberra
treatment_group	1.000	0.001	0.001	0.608	≥0.25	ns
Residual	71.000	0.150	0.002	NA	NA	NA

Tukey: Homogeneity of Dispersion
Tukey(Beta Disperson ~ Treatment); Parasite exposure only 60 DPE
.y.	term	group1	group2	estimate	conf.low	conf.high	adj.p.value	sig
bray
Distance	treatment_group	A- T- P+	A- T- P-	−0.025	−0.105	0.055	≥0.25	ns
canberra
Distance	treatment_group	A- T- P+	A- T- P-	0.008	−0.013	0.030	≥0.25	ns

CAP

Plot

Code

# Set the subsection for this analysis
tmp.resSubSection <- "ParasiteExp_60DPE"


#### Capscale ----------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.mod"]] <-
  run_capscale(ps.cleaned_alpha %>% # is the phyloseq object we are looping over
          microViz::ps_filter(Treatment  %in% c(
          "A- T- P-",
          "A- T- P+"#,
          # "A+ T- P-",
          # "A+ T- P+",
          # "A- T+ P-",
          # "A- T+ P+",
          # "A+ T+ P-",
          # "A+ T+ P+"
          )) %>%
          microViz::ps_filter(Time == 60), 
               dist.matrix = beta.dist.mat[[tmp.resSubSection]], 
               formula_str = "dist ~ Treatment")

##### ADONIS ------------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS"]] <-
  run_cap_adonis(ps.cleaned_alpha %>% # is the phyloseq object we are looping over
          microViz::ps_filter(Treatment  %in% c(
          "A- T- P-",
          "A- T- P+"#,
          # "A+ T- P-",
          # "A+ T- P+",
          # "A- T+ P-",
          # "A- T+ P+",
          # "A+ T+ P-",
          # "A+ T+ P+"
          )) %>%
          microViz::ps_filter(Time == 60),
                 dist.matrix = beta.dist.mat[[tmp.resSubSection]], 
                 formula_str = "dist ~ Treatment",
                 by.method = "terms") 

## Warning in tidy.anova(.): The following column names in ANOVA output were not
## recognized or transformed: SumOfSqs, R2

## Warning in tidy.anova(.): The following column names in ANOVA output were not
## recognized or transformed: SumOfSqs, R2

### Table

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS.Table"]] <-
  beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS"]] %>%
    set_GT(var = "p.value", group.by = "Beta.Metric") %>%
  
  # Title/caption
  gt::tab_header(
    title = "ADONIS2",
    subtitle = "adonis2(Beta Distance ~ Treatment); Parasite exposure only at 60 DPE"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = c("p.value", "sig"),
      rows = p.value < 0.05
    )
  )

Tables

ADONIS2
adonis2(Beta Distance ~ Treatment); Parasite exposure only at 60 DPE
term	df	SumOfSqs	R2	statistic	p.value	sig
bray
Treatment	1.000	0.587	0.047	3.468	0.007	**
Residual	71.000	12.011	0.953	NA	NA	NA
Total	72.000	12.597	1.000	NA	NA	NA
canberra
Treatment	1.000	1.287	0.050	3.741	0.001	***
Residual	71.000	24.425	0.950	NA	NA	NA
Total	72.000	25.712	1.000	NA	NA	NA

02) Is the host-microbiome response to parasite exposure historically contingent?

Alpha

All

Displacement

Plot

Code

# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressParaExp"

# GLM Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]] <- 
  test.alpha.all_metrics$volatility %>%
    # Filter 
    dplyr::filter(Treatment  %in% c(
          # "A- T- P-",
          "A- T- P+",
          # "A+ T- P-",
          "A+ T- P+",
          # "A- T+ P-",
          "A- T+ P+",
          # "A+ T+ P-",
          "A+ T+ P+"
          )) %>%
    
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models(formula_str = "Displacement ~ Treatment", 
                   alpha_metric_col = "Alpha.Metric", 
                   alpha_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["ANOVA"]] <-
  run_glm_anova(alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]])

# Tukey Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["Tukey"]] <-
  test.alpha.all_metrics$volatility %>%
    # Filter 
    dplyr::filter(Treatment  %in% c(
          # "A- T- P-",
          "A- T- P+",
          # "A+ T- P-",
          "A+ T- P+",
          # "A- T+ P-",
          "A- T+ P+",
          # "A+ T+ P-",
          "A+ T+ P+"
          )) %>%
    
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm(., "Displacement", "Alpha.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 

Tables

GLM Results
glm(Displacement ~ Treatment); Prior stressors and parasite exposure
term	estimate	std.error	statistic	p.value	p.adj.sig
Shannon
(Intercept)	0.033	0.020	1.674	0.095	ns
A+ T- P+	−0.027	0.028	−0.956	≥0.25	ns
A- T+ P+	−0.047	0.028	−1.675	0.095	ns
A+ T+ P+	−0.017	0.029	−0.575	≥0.25	ns
Simpson
(Intercept)	−0.006	0.028	−0.219	≥0.25	ns
A+ T- P+	−0.037	0.040	−0.916	≥0.25	ns
A- T+ P+	−0.105	0.041	−2.592	0.010	*
A+ T+ P+	−0.042	0.042	−1.020	≥0.25	ns
Richness
(Intercept)	0.230	0.025	9.206	<0.001	****
A+ T- P+	−0.054	0.036	−1.519	0.130	ns
A- T+ P+	0.007	0.036	0.186	≥0.25	ns
A+ T+ P+	−0.063	0.037	−1.711	0.088	ns

ANOVA of GLM
ANOVA(GLM(Displacement ~ Treatment), type = 2); Prior stressors and parasite exposure
term	statistic	df	p.value	sig
Shannon
Treatment	2.931	3.000	≥0.25	ns
Simpson
Treatment	6.896	3.000	0.075	ns
Richness
Treatment	5.804	3.000	0.122	ns

Pairwise Tukey's HSD, p.adj: Dunnett
Tukey(Displacement ~ Treatment); Prior stressors and parasite exposure
term	.y.	group1	group2	estimate	std.error	statistic	adj.p.value	Variable
Shannon
Treatment	Alpha.Score	A+ T- P+	A- T- P+	−0.027	0.028	−0.956	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A- T- P+	−0.047	0.028	−1.675	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T- P+	−0.017	0.029	−0.575	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A+ T- P+	−0.020	0.028	−0.721	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A+ T- P+	0.010	0.029	0.349	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T+ P+	0.031	0.029	1.046	≥0.25	Treatment
Simpson
Treatment	Alpha.Score	A+ T- P+	A- T- P+	−0.037	0.040	−0.916	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A- T- P+	−0.105	0.041	−2.592	0.047	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T- P+	−0.042	0.042	−1.020	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A+ T- P+	−0.068	0.041	−1.671	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A+ T- P+	−0.005	0.042	−0.131	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T+ P+	0.063	0.042	1.490	≥0.25	Treatment
Richness
Treatment	Alpha.Score	A+ T- P+	A- T- P+	−0.054	0.036	−1.519	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A- T- P+	0.007	0.036	0.186	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T- P+	−0.063	0.037	−1.711	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A+ T- P+	0.061	0.036	1.685	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A+ T- P+	−0.009	0.037	−0.237	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T+ P+	−0.070	0.037	−1.870	0.241	Treatment

60 DPE

Displacement

Plot

Code

# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressParaExp_60DPE"

# GLM Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]] <- 
  test.alpha.all_metrics$displacement %>%
    dplyr::filter(Treatment  %in% c(
          # "A- T- P-",
          "A- T- P+",
          # "A+ T- P-",
          "A+ T- P+",
          # "A- T+ P-",
          "A- T+ P+",
          # "A+ T+ P-",
          "A+ T+ P+"
          )) %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models(formula_str = "Displacement ~ Treatment", 
                   alpha_metric_col = "Alpha.Metric", 
                   alpha_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["ANOVA"]] <-
  run_glm_anova(alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]])

# Tukey Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["Tukey"]] <-
  test.alpha.all_metrics$displacement %>%
    dplyr::filter(Treatment  %in% c(
          # "A- T- P-",
          "A- T- P+",
          # "A+ T- P-",
          "A+ T- P+",
          # "A- T+ P-",
          "A- T+ P+",
          # "A+ T+ P-",
          "A+ T+ P+"
          )) %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm(., "Displacement", "Alpha.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 

Tables

GLM Results - 60 DPE
glm(Displacement ~ Treatment); All treatments at 60 DPE
term	estimate	std.error	statistic	p.value	p.adj.sig
Shannon
(Intercept)	0.008	0.037	0.224	≥0.25	ns
A+ T- P+	−0.001	0.053	−0.022	≥0.25	ns
A- T+ P+	−0.077	0.054	−1.424	0.157	ns
A+ T+ P+	0.004	0.058	0.065	≥0.25	ns
Simpson
(Intercept)	−0.062	0.048	−1.286	0.201	ns
A+ T- P+	0.000	0.069	−0.005	≥0.25	ns
A- T+ P+	−0.180	0.071	−2.547	0.012	*
A+ T+ P+	−0.037	0.076	−0.484	≥0.25	ns
Richness
(Intercept)	0.363	0.032	11.248	<0.001	****
A+ T- P+	−0.132	0.046	−2.863	0.005	**
A- T+ P+	−0.003	0.047	−0.069	≥0.25	ns
A+ T+ P+	−0.053	0.051	−1.039	≥0.25	ns

ANOVA of GLM - 60 DPE
ANOVA(GLM(Displacement ~ Treatment), type = 2); Prior stressors and parasite exposure at 60 DPE
term	statistic	df	p.value	sig
Shannon
Treatment	2.892	3.000	≥0.25	ns
Simpson
Treatment	8.445	3.000	0.038	*
Richness
Treatment	10.375	3.000	0.016	*

Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); Prior stressors and parasite exposure at 60 DPE
term	.y.	group1	group2	estimate	std.error	statistic	adj.p.value	Variable
Shannon
Treatment	Alpha.Score	A+ T- P+	A- T- P+	−0.001	0.053	−0.022	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A- T- P+	−0.077	0.054	−1.424	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T- P+	0.004	0.058	0.065	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A+ T- P+	−0.076	0.055	−1.384	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A+ T- P+	0.005	0.059	0.084	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T+ P+	0.081	0.060	1.351	≥0.25	Treatment
Simpson
Treatment	Alpha.Score	A+ T- P+	A- T- P+	0.000	0.069	−0.005	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A- T- P+	−0.180	0.071	−2.547	0.053	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T- P+	−0.037	0.076	−0.484	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A+ T- P+	−0.179	0.071	−2.509	0.058	Treatment
Treatment	Alpha.Score	A+ T+ P+	A+ T- P+	−0.036	0.077	−0.475	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T+ P+	0.143	0.078	1.834	≥0.25	Treatment
Richness
Treatment	Alpha.Score	A+ T- P+	A- T- P+	−0.132	0.046	−2.863	0.022	Treatment
Treatment	Alpha.Score	A- T+ P+	A- T- P+	−0.003	0.047	−0.069	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T- P+	−0.053	0.051	−1.039	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P+	A+ T- P+	0.129	0.048	2.692	0.036	Treatment
Treatment	Alpha.Score	A+ T+ P+	A+ T- P+	0.080	0.051	1.551	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P+	A- T+ P+	−0.050	0.052	−0.946	≥0.25	Treatment

Significant Results:Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); All treatments at 60 DPE
term	.y.	group1	group2	estimate	std.error	statistic	adj.p.value	Variable
Richness
Treatment	Alpha.Score	A+ T- P+	A- T- P+	−0.132	0.046	−2.863	0.022	Treatment
Treatment	Alpha.Score	A- T+ P+	A+ T- P+	0.129	0.048	2.692	0.036	Treatment

Beta

All

Displacement

Plot

Code

# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressParaExp"

# GLM Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]] <- 
  test.beta.all_metrics$volatility %>%
        # Filter 
    dplyr::filter(Treatment  %in% c(
          # "A- T- P-",
          "A- T- P+",
          # "A+ T- P-",
          "A+ T- P+",
          # "A- T+ P-",
          "A- T+ P+",
          # "A+ T+ P-",
          "A+ T+ P+"
          )) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models_beta(formula_str = "Displacement ~ Treatment", 
                   beta_metric_col = "Beta.Metric", 
                   beta_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["ANOVA"]] <-
  run_glm_anova_beta(beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]])

# Tukey Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["Tukey"]] <-
  test.beta.all_metrics$volatility %>%
        # Filter 
    dplyr::filter(Treatment  %in% c(
          # "A- T- P-",
          "A- T- P+",
          # "A+ T- P-",
          "A+ T- P+",
          # "A- T+ P-",
          "A- T+ P+",
          # "A+ T+ P-",
          "A+ T+ P+"
          )) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm_beta(., "Displacement", "Beta.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 

Tables

GLM Results - Beta
glm(Displacement ~ Treatment); Prior stressors and parasite exposure
term	estimate	std.error	statistic	p.value	p.adj.sig
Bray-Curtis_norm
(Intercept)	0.081	0.027	2.973	0.003	**
A+ T- P+	−0.141	0.039	−3.633	<0.001	***
A- T+ P+	−0.306	0.039	−7.843	<0.001	****
A+ T+ P+	−0.146	0.040	−3.662	<0.001	***
Canberra_norm
(Intercept)	−0.106	0.020	−5.372	<0.001	****
A+ T- P+	0.013	0.028	0.478	≥0.25	ns
A- T+ P+	−0.016	0.028	−0.573	≥0.25	ns
A+ T+ P+	0.026	0.029	0.900	≥0.25	ns

ANOVA of GLM - Beta
ANOVA(GLM(Displacement ~ Treatment), type = 2); Prior stressors and parasite exposure
term	statistic	df	p.value	sig
Bray-Curtis_norm
Treatment	61.580	3.000	<0.001	****
Canberra_norm
Treatment	2.326	3.000	≥0.25	ns

Pairwise Tukey's HSD, p.adj: Dunnett - Beta
Tukey(Displacement ~ Treatment); Prior stressors and parasite exposure
term	.y.	group1	group2	estimate	std.error	statistic	adj.p.value	Variable
Bray-Curtis_norm
Treatment	Beta.Score	A+ T- P+	A- T- P+	−0.141	0.039	−3.633	0.002	Treatment
Treatment	Beta.Score	A- T+ P+	A- T- P+	−0.306	0.039	−7.843	<0.001	Treatment
Treatment	Beta.Score	A+ T+ P+	A- T- P+	−0.146	0.040	−3.662	0.002	Treatment
Treatment	Beta.Score	A- T+ P+	A+ T- P+	−0.165	0.039	−4.208	<0.001	Treatment
Treatment	Beta.Score	A+ T+ P+	A+ T- P+	−0.006	0.040	−0.141	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P+	A- T+ P+	0.159	0.040	3.940	<0.001	Treatment
Canberra_norm
Treatment	Beta.Score	A+ T- P+	A- T- P+	0.013	0.028	0.478	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P+	A- T- P+	−0.016	0.028	−0.573	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P+	A- T- P+	0.026	0.029	0.900	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P+	A+ T- P+	−0.030	0.029	−1.042	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P+	A+ T- P+	0.013	0.029	0.434	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P+	A- T+ P+	0.042	0.029	1.442	≥0.25	Treatment

60 DPE

Displacement

Plot

Code

# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressParaExp_60DPE"

# GLM Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]] <- 
  test.beta.all_metrics$displacement %>%
    dplyr::filter(Treatment  %in% c(
          # "A- T- P-",
          "A- T- P+",
          # "A+ T- P-",
          "A+ T- P+",
          # "A- T+ P-",
          "A- T+ P+",
          # "A+ T+ P-",
          "A+ T+ P+"
          )) %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models_beta(formula_str = "Displacement ~ Treatment", 
                   beta_metric_col = "Beta.Metric", 
                   beta_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["ANOVA"]] <-
  run_glm_anova_beta(beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]])

# Tukey Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["Tukey"]] <-
  test.beta.all_metrics$displacement %>%
    dplyr::filter(Treatment  %in% c(
          # "A- T- P-",
          "A- T- P+",
          # "A+ T- P-",
          "A+ T- P+",
          # "A- T+ P-",
          "A- T+ P+",
          # "A+ T+ P-",
          "A+ T+ P+"
          )) %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm_beta(., "Displacement", "Beta.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 

Tables

GLM Results - 60 DPE
glm(Displacement ~ Treatment); Prior stressors and parasite exposure only at 60 DPE
term	estimate	std.error	statistic	p.value	p.adj.sig
Bray-Curtis_norm
(Intercept)	0.224	0.038	5.823	<0.001	****
A+ T- P+	−0.305	0.055	−5.534	<0.001	****
A- T+ P+	−0.490	0.056	−8.688	<0.001	****
A+ T+ P+	−0.220	0.060	−3.646	<0.001	***
Canberra_norm
(Intercept)	0.010	0.028	0.340	≥0.25	ns
A+ T- P+	−0.097	0.040	−2.398	0.018	*
A- T+ P+	−0.069	0.041	−1.670	0.097	ns
A+ T+ P+	−0.063	0.044	−1.425	0.157	ns

ANOVA of GLM - 60 DPE
ANOVA(GLM(Displacement ~ Treatment), type = 2); Prior stressors and parasite exposure only at 60 DPE
term	statistic	df	p.value	sig
Bray-Curtis_norm
Treatment	78.442	3.000	<0.001	****
Canberra_norm
Treatment	6.174	3.000	0.103	ns

Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); Prior stressors and parasite exposure only at 60 DPE
term	.y.	group1	group2	estimate	std.error	statistic	adj.p.value	Variable
Bray-Curtis_norm
Treatment	Beta.Score	A+ T- P+	A- T- P+	−0.305	0.055	−5.534	<0.001	Treatment
Treatment	Beta.Score	A- T+ P+	A- T- P+	−0.490	0.056	−8.688	0.000	Treatment
Treatment	Beta.Score	A+ T+ P+	A- T- P+	−0.220	0.060	−3.646	0.002	Treatment
Treatment	Beta.Score	A- T+ P+	A+ T- P+	−0.185	0.057	−3.240	0.007	Treatment
Treatment	Beta.Score	A+ T+ P+	A+ T- P+	0.084	0.061	1.379	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P+	A- T+ P+	0.269	0.062	4.321	<0.001	Treatment
Canberra_norm
Treatment	Beta.Score	A+ T- P+	A- T- P+	−0.097	0.040	−2.398	0.077	Treatment
Treatment	Beta.Score	A- T+ P+	A- T- P+	−0.069	0.041	−1.670	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P+	A- T- P+	−0.063	0.044	−1.425	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P+	A+ T- P+	0.028	0.042	0.663	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P+	A+ T- P+	0.034	0.045	0.750	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P+	A- T+ P+	0.006	0.046	0.128	≥0.25	Treatment

Significant Results:Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); All treatments at 60 DPE
term	.y.	group1	group2	estimate	std.error	statistic	adj.p.value	Variable
Bray-Curtis_norm
Treatment	Beta.Score	A+ T- P+	A- T- P+	−0.305	0.055	−5.534	<0.001	Treatment
Treatment	Beta.Score	A- T+ P+	A- T- P+	−0.490	0.056	−8.688	0.000	Treatment
Treatment	Beta.Score	A+ T+ P+	A- T- P+	−0.220	0.060	−3.646	0.002	Treatment
Treatment	Beta.Score	A- T+ P+	A+ T- P+	−0.185	0.057	−3.240	0.007	Treatment
Treatment	Beta.Score	A+ T+ P+	A- T+ P+	0.269	0.062	4.321	<0.001	Treatment

PCoA

Plot

Code

# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressParaExp_60DPE"



#### Dispersion --------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]] <-
  run_BetaDispersion(dist.matrix = beta.dist.mat[[tmp.resSubSection]], 
                     var = c("treatment_group"))

##### ANOVA ------------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA"]] <-
  get_HoD_anova(betaDisper = beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]],
                var = c("treatment_group"))

### Table

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA.Table"]] <-
  beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA"]] %>%
  
    # Create the table
    dplyr::group_by(Beta.Metric) %>%
    set_GT(var = "p.value", group.by = "Beta.Metric")  %>%
  
    # Title/caption
    gt::tab_header(
      title = "ANOVA: Homogeneity of Dispersion",
      subtitle = "ANOVA(Beta Disperson ~ Treatment); Prior stressors and parasite exposure at 60 DPE"
    ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = c("p.value", "sig"),
      rows = p.value < 0.05
    )
  )

##### Tukey ------------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey"]] <-
  get_HoD_tukey(betaDisper = beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]],
                var = c("treatment_group")) 

### Table

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey.Table"]] <-
  beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey"]] %>%
    
    dplyr::mutate(
        group1 = case_when(
        group1 == "Control" ~ "A- T- P-",
        group1 == "Parasite" ~ "A- T- P+",
        group1 == "Antibiotics" ~ "A+ T- P-",
        group1 == "Antibiotics_Parasite" ~ "A+ T- P+",
        group1 == "Temperature" ~ "A- T+ P-",
        group1 == "Temperature_Parasite" ~ "A- T+ P+",
        group1 == "Antibiotics_Temperature" ~ "A+ T+ P-",
        group1 == "Antibiotics_Temperature_Parasite" ~ "A+ T+ P+",
    ),
        group2 = case_when(
        group2 == "Control" ~ "A- T- P-",
        group2 == "Parasite" ~ "A- T- P+",
        group2 == "Antibiotics" ~ "A+ T- P-",
        group2 == "Antibiotics_Parasite" ~ "A+ T- P+",
        group2 == "Temperature" ~ "A- T+ P-",
        group2 == "Temperature_Parasite" ~ "A- T+ P+",
        group2 == "Antibiotics_Temperature" ~ "A+ T+ P-",
        group2 == "Antibiotics_Temperature_Parasite" ~ "A+ T+ P+",
    )) %>%
  
  # Create the table
  dplyr::group_by(Beta.Metric) %>%
  set_GT(var = "adj.p.value", group.by = "Beta.Metric")  %>%
  
  # Title/caption
  gt::tab_header(
    title = "Tukey: Homogeneity of Dispersion",
    subtitle = "Tukey(Beta Disperson ~ Treatment); Prior stressors and parasite exposure at 60 DPE"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = c("adj.p.value", "sig"),
      rows = adj.p.value < 0.05
    )
  )

Tables

ANOVA: Homogeneity of Dispersion
ANOVA(Beta Disperson ~ Treatment); Prior stressors and parasite exposure at 60 DPE
term	df	sumsq	meansq	statistic	p.value	sig
bray
treatment_group	3.000	0.291	0.097	3.785	0.012	*
Residual	125.000	3.198	0.026	NA	NA	NA
canberra
treatment_group	3.000	0.049	0.016	4.464	0.005	**
Residual	125.000	0.462	0.004	NA	NA	NA

Tukey: Homogeneity of Dispersion
Tukey(Beta Disperson ~ Treatment); Prior stressors and parasite exposure at 60 DPE
.y.	term	group1	group2	estimate	conf.low	conf.high	adj.p.value	sig
bray
Distance	treatment_group	A+ T- P+	A- T- P+	−0.040	−0.138	0.058	≥0.25	ns
Distance	treatment_group	A- T+ P+	A- T- P+	−0.081	−0.182	0.019	0.157	ns
Distance	treatment_group	A+ T+ P+	A- T- P+	−0.132	−0.240	−0.024	0.009	**
Distance	treatment_group	A- T+ P+	A+ T- P+	−0.041	−0.143	0.061	≥0.25	ns
Distance	treatment_group	A+ T+ P+	A+ T- P+	−0.092	−0.201	0.017	0.130	ns
Distance	treatment_group	A+ T+ P+	A- T+ P+	−0.051	−0.162	0.060	≥0.25	ns
canberra
Distance	treatment_group	A+ T- P+	A- T- P+	−0.001	−0.038	0.036	≥0.25	ns
Distance	treatment_group	A- T+ P+	A- T- P+	−0.015	−0.053	0.023	≥0.25	ns
Distance	treatment_group	A+ T+ P+	A- T- P+	−0.052	−0.093	−0.011	0.007	**
Distance	treatment_group	A- T+ P+	A+ T- P+	−0.014	−0.053	0.025	≥0.25	ns
Distance	treatment_group	A+ T+ P+	A+ T- P+	−0.051	−0.093	−0.010	0.009	**
Distance	treatment_group	A+ T+ P+	A- T+ P+	−0.037	−0.079	0.005	0.106	ns

CAP

Plot

Code

# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressParaExp_60DPE"


#### Capscale ----------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.mod"]] <-
  run_capscale(ps.cleaned_alpha %>% # is the phyloseq object we are looping over
          microViz::ps_filter(Treatment  %in% c(
          # "A- T- P-",
          "A- T- P+",
          # "A+ T- P-",
          "A+ T- P+",
          # "A- T+ P-",
          "A- T+ P+",
          # "A+ T+ P-",
          "A+ T+ P+"
          )) %>%
          microViz::ps_filter(Time == 60), 
               dist.matrix = beta.dist.mat[[tmp.resSubSection]], 
               formula_str = "dist ~ Treatment")

##### ADONIS ------------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS"]] <-
  run_cap_adonis(ps.cleaned_alpha %>% # is the phyloseq object we are looping over
          microViz::ps_filter(Treatment  %in% c(
          # "A- T- P-",
          "A- T- P+",
          # "A+ T- P-",
          "A+ T- P+",
          # "A- T+ P-",
          "A- T+ P+",
          # "A+ T+ P-",
          "A+ T+ P+"
          )) %>%
          microViz::ps_filter(Time == 60),
                 dist.matrix = beta.dist.mat[[tmp.resSubSection]], 
                 formula_str = "dist ~ Treatment",
                 by.method = "terms") 

## Warning in tidy.anova(.): The following column names in ANOVA output were not
## recognized or transformed: SumOfSqs, R2

## Warning in tidy.anova(.): The following column names in ANOVA output were not
## recognized or transformed: SumOfSqs, R2

### Table

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS.Table"]] <-
  beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS"]] %>%
    set_GT(var = "p.value", group.by = "Beta.Metric") %>%
  
  # Title/caption
  gt::tab_header(
    title = "ADONIS2",
    subtitle = "adonis2(Beta Distance ~ Treatment); Prior stressors and parasite exposure at 60 DPE"
  )  %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = c("p.value", "sig"),
      rows = p.value < 0.05
    )
  )

Tables

ADONIS2
adonis2(Beta Distance ~ Treatment); Prior stressors and parasite exposure at 60 DPE
term	df	SumOfSqs	R2	statistic	p.value	sig
bray
Treatment	3.000	0.708	0.045	1.959	0.037	*
Residual	125.000	15.053	0.955	NA	NA	NA
Total	128.000	15.761	1.000	NA	NA	NA
canberra
Treatment	3.000	2.852	0.064	2.833	0.001	***
Residual	125.000	41.942	0.936	NA	NA	NA
Total	128.000	44.794	1.000	NA	NA	NA

03) Is the host-microbiome recovery historically contingent?

Alpha

All

Displacement

Plot

Code

# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressNoParaExp"

# GLM Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]] <- 
  test.alpha.all_metrics$volatility %>%
    # Filter 
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          # "A- T- P+",
          "A+ T- P-",
          # "A+ T- P+",
          "A- T+ P-",
          # "A- T+ P+",
          "A+ T+ P-"#,
          # "A+ T+ P+"
          )) %>%
    
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models(formula_str = "Displacement ~ Treatment", 
                   alpha_metric_col = "Alpha.Metric", 
                   alpha_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["ANOVA"]] <-
  run_glm_anova(alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]])

# Tukey Analysis
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["Tukey"]] <-
  test.alpha.all_metrics$volatility %>%
    # Filter 
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          # "A- T- P+",
          "A+ T- P-",
          # "A+ T- P+",
          "A- T+ P-",
          # "A- T+ P+",
          "A+ T+ P-"#,
          # "A+ T+ P+"
          )) %>%
    
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm(., "Displacement", "Alpha.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 

Tables

GLM Results
glm(Displacement ~ Treatment); Prior stressors and no parasite exposure
term	estimate	std.error	statistic	p.value	p.adj.sig
Shannon
(Intercept)	0.078	0.019	4.047	<0.001	****
A+ T- P-	−0.051	0.028	−1.839	0.067	ns
A- T+ P-	−0.004	0.028	−0.153	≥0.25	ns
A+ T+ P-	−0.113	0.029	−3.849	<0.001	***
Simpson
(Intercept)	0.042	0.027	1.562	0.119	ns
A+ T- P-	−0.034	0.039	−0.861	≥0.25	ns
A- T+ P-	0.035	0.040	0.894	≥0.25	ns
A+ T+ P-	−0.181	0.042	−4.369	<0.001	****
Richness
(Intercept)	0.221	0.023	9.748	<0.001	****
A+ T- P-	−0.049	0.033	−1.508	0.133	ns
A- T+ P-	−0.016	0.033	−0.493	≥0.25	ns
A+ T+ P-	0.029	0.035	0.825	≥0.25	ns

ANOVA of GLM
ANOVA(GLM(Displacement ~ Treatment), type = 2); Prior stressors and no parasite exposure
term	statistic	df	p.value	sig
Shannon
Treatment	18.423	3.000	<0.001	***
Simpson
Treatment	29.201	3.000	<0.001	****
Richness
Treatment	5.193	3.000	0.158	ns

Pairwise Tukey's HSD, p.adj: Dunnett
Tukey(Displacement ~ Treatment); Prior stressors and no parasite exposure
term	.y.	group1	group2	estimate	std.error	statistic	adj.p.value	Variable
Shannon
Treatment	Alpha.Score	A+ T- P-	A- T- P-	−0.051	0.028	−1.839	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P-	A- T- P-	−0.004	0.028	−0.153	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T- P-	−0.113	0.029	−3.849	<0.001	Treatment
Treatment	Alpha.Score	A- T+ P-	A+ T- P-	0.047	0.029	1.629	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A+ T- P-	−0.062	0.030	−2.083	0.158	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T+ P-	−0.109	0.030	−3.596	0.002	Treatment
Simpson
Treatment	Alpha.Score	A+ T- P-	A- T- P-	−0.034	0.039	−0.861	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P-	A- T- P-	0.035	0.040	0.894	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T- P-	−0.181	0.042	−4.369	<0.001	Treatment
Treatment	Alpha.Score	A- T+ P-	A+ T- P-	0.069	0.040	1.711	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A+ T- P-	−0.148	0.042	−3.500	0.003	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T+ P-	−0.217	0.043	−5.073	<0.001	Treatment
Richness
Treatment	Alpha.Score	A+ T- P-	A- T- P-	−0.049	0.033	−1.508	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P-	A- T- P-	−0.016	0.033	−0.493	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T- P-	0.029	0.035	0.825	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P-	A+ T- P-	0.033	0.034	0.975	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A+ T- P-	0.078	0.035	2.206	0.121	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T+ P-	0.045	0.036	1.258	≥0.25	Treatment

60 DPE

Displacement

Plot

Code

# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressNoParaExp_60DPE"

# GLM Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]] <- 
  test.alpha.all_metrics$displacement %>%
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          # "A- T- P+",
          "A+ T- P-",
          # "A+ T- P+",
          "A- T+ P-",
          # "A- T+ P+",
          "A+ T+ P-"#,
          # "A+ T+ P+"
          )) %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models(formula_str = "Displacement ~ Treatment", 
                   alpha_metric_col = "Alpha.Metric", 
                   alpha_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["ANOVA"]] <-
  run_glm_anova(alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]])

# Tukey Analysis for 60 DPE
alpha.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["Tukey"]] <-
  test.alpha.all_metrics$displacement %>%
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          # "A- T- P+",
          "A+ T- P-",
          # "A+ T- P+",
          "A- T+ P-",
          # "A- T+ P+",
          "A+ T+ P-"#,
          # "A+ T+ P+"
          )) %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Alpha.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm(., "Displacement", "Alpha.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 

Tables

GLM Results - 60 DPE
glm(Displacement ~ Treatment); Prior stressors and no parasite exposure at 60 DPE
term	estimate	std.error	statistic	p.value	p.adj.sig
Shannon
(Intercept)	0.091	0.036	2.496	0.014	*
A+ T- P-	−0.034	0.054	−0.627	≥0.25	ns
A- T+ P-	−0.041	0.056	−0.731	≥0.25	ns
A+ T+ P-	−0.191	0.067	−2.860	0.005	**
Simpson
(Intercept)	0.008	0.050	0.161	≥0.25	ns
A+ T- P-	0.023	0.074	0.308	≥0.25	ns
A- T+ P-	0.030	0.077	0.395	≥0.25	ns
A+ T+ P-	−0.332	0.092	−3.627	<0.001	***
Richness
(Intercept)	0.338	0.031	10.867	<0.001	****
A+ T- P-	−0.064	0.046	−1.388	0.168	ns
A- T+ P-	−0.089	0.048	−1.856	0.066	ns
A+ T+ P-	0.067	0.057	1.174	0.243	ns

ANOVA of GLM - 60 DPE
ANOVA(GLM(Displacement ~ Treatment), type = 2); Prior stressors and no parasite exposure at 60 DPE
term	statistic	df	p.value	sig
Shannon
Treatment	8.403	3.000	0.038	*
Simpson
Treatment	17.775	3.000	<0.001	***
Richness
Treatment	8.617	3.000	0.035	*

Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); Prior stressors and no parasite exposure at 60 DPE
term	.y.	group1	group2	estimate	std.error	statistic	adj.p.value	Variable
Shannon
Treatment	Alpha.Score	A+ T- P-	A- T- P-	−0.034	0.054	−0.627	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P-	A- T- P-	−0.041	0.056	−0.731	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T- P-	−0.191	0.067	−2.860	0.022	Treatment
Treatment	Alpha.Score	A- T+ P-	A+ T- P-	−0.007	0.058	−0.124	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A+ T- P-	−0.158	0.069	−2.289	0.099	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T+ P-	−0.150	0.071	−2.131	0.142	Treatment
Simpson
Treatment	Alpha.Score	A+ T- P-	A- T- P-	0.023	0.074	0.308	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P-	A- T- P-	0.030	0.077	0.395	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T- P-	−0.332	0.092	−3.627	0.002	Treatment
Treatment	Alpha.Score	A- T+ P-	A+ T- P-	0.008	0.080	0.095	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A+ T- P-	−0.355	0.094	−3.765	<0.001	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T+ P-	−0.362	0.097	−3.751	0.001	Treatment
Richness
Treatment	Alpha.Score	A+ T- P-	A- T- P-	−0.064	0.046	−1.388	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P-	A- T- P-	−0.089	0.048	−1.856	0.245	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T- P-	0.067	0.057	1.174	≥0.25	Treatment
Treatment	Alpha.Score	A- T+ P-	A+ T- P-	−0.025	0.050	−0.502	≥0.25	Treatment
Treatment	Alpha.Score	A+ T+ P-	A+ T- P-	0.131	0.059	2.226	0.115	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T+ P-	0.156	0.061	2.586	0.047	Treatment

Significant Results:Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); All treatments at 60 DPE
term	.y.	group1	group2	estimate	std.error	statistic	adj.p.value	Variable
Shannon
Treatment	Alpha.Score	A+ T+ P-	A- T- P-	−0.191	0.067	−2.860	0.022	Treatment
Simpson
Treatment	Alpha.Score	A+ T+ P-	A- T- P-	−0.332	0.092	−3.627	0.002	Treatment
Treatment	Alpha.Score	A+ T+ P-	A+ T- P-	−0.355	0.094	−3.765	<0.001	Treatment
Treatment	Alpha.Score	A+ T+ P-	A- T+ P-	−0.362	0.097	−3.751	0.001	Treatment
Richness
Treatment	Alpha.Score	A+ T+ P-	A- T+ P-	0.156	0.061	2.586	0.047	Treatment

Beta

All

Displacement

Plot

Code

# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressNoParaExp"

# GLM Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]] <- 
  test.beta.all_metrics$volatility %>%
        # Filter 
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          # "A- T- P+",
          "A+ T- P-",
          # "A+ T- P+",
          "A- T+ P-",
          # "A- T+ P+",
          "A+ T+ P-"#,
          # "A+ T+ P+"
          )) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models_beta(formula_str = "Displacement ~ Treatment", 
                   beta_metric_col = "Beta.Metric", 
                   beta_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["ANOVA"]] <-
  run_glm_anova_beta(beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["GLM"]])

# Tukey Analysis for Beta
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["ALL"]][["Tukey"]] <-
  test.beta.all_metrics$volatility %>%
        # Filter 
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          # "A- T- P+",
          "A+ T- P-",
          # "A+ T- P+",
          "A- T+ P-",
          # "A- T+ P+",
          "A+ T+ P-"#,
          # "A+ T+ P+"
          )) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm_beta(., "Displacement", "Beta.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 

Tables

GLM Results - Beta
glm(Displacement ~ Treatment); Prior stressors and no parasite exposure
term	estimate	std.error	statistic	p.value	p.adj.sig
Bray-Curtis_norm
(Intercept)	−0.061	0.027	−2.279	0.023	*
A+ T- P-	0.016	0.039	0.414	≥0.25	ns
A- T+ P-	0.081	0.039	2.059	0.040	*
A+ T+ P-	−0.020	0.041	−0.475	≥0.25	ns
Canberra_norm
(Intercept)	−0.176	0.023	−7.610	<0.001	****
A+ T- P-	0.040	0.033	1.213	0.226	ns
A- T+ P-	0.179	0.034	5.285	<0.001	****
A+ T+ P-	0.155	0.035	4.359	<0.001	****

ANOVA of GLM - Beta
ANOVA(GLM(Displacement ~ Treatment), type = 2); Prior stressors and no parasite exposure
term	statistic	df	p.value	sig
Bray-Curtis_norm
Treatment	6.694	3.000	0.082	ns
Canberra_norm
Treatment	38.039	3.000	<0.001	****

Pairwise Tukey's HSD, p.adj: Dunnett - Beta
Tukey(Displacement ~ Treatment); Prior stressors and no parasite exposure
term	.y.	group1	group2	estimate	std.error	statistic	adj.p.value	Variable
Bray-Curtis_norm
Treatment	Beta.Score	A+ T- P-	A- T- P-	0.016	0.039	0.414	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P-	A- T- P-	0.081	0.039	2.059	0.166	Treatment
Treatment	Beta.Score	A+ T+ P-	A- T- P-	−0.020	0.041	−0.475	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P-	A+ T- P-	0.065	0.040	1.621	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P-	A+ T- P-	−0.036	0.042	−0.851	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P-	A- T+ P-	−0.101	0.042	−2.371	0.083	Treatment
Canberra_norm
Treatment	Beta.Score	A+ T- P-	A- T- P-	0.040	0.033	1.213	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P-	A- T- P-	0.179	0.034	5.285	<0.001	Treatment
Treatment	Beta.Score	A+ T+ P-	A- T- P-	0.155	0.035	4.359	<0.001	Treatment
Treatment	Beta.Score	A- T+ P-	A+ T- P-	0.138	0.034	4.015	<0.001	Treatment
Treatment	Beta.Score	A+ T+ P-	A+ T- P-	0.114	0.036	3.164	0.009	Treatment
Treatment	Beta.Score	A+ T+ P-	A- T+ P-	−0.024	0.037	−0.666	≥0.25	Treatment

60 DPE

Displacement

Plot

Code

# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressNoParaExp_60DPE"

# GLM Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]] <- 
  test.beta.all_metrics$displacement %>%
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          # "A- T- P+",
          "A+ T- P-",
          # "A+ T- P+",
          "A- T+ P-",
          # "A- T+ P+",
          "A+ T+ P-"#,
          # "A+ T+ P+"
          )) %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run GLM models
  run_glm_models_beta(formula_str = "Displacement ~ Treatment", 
                   beta_metric_col = "Beta.Metric", 
                   beta_score_col = "Displacement",
                   family_str = "gaussian")

# ANOVA Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["ANOVA"]] <-
  run_glm_anova_beta(beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["GLM"]])

# Tukey Analysis for 60 DPE
beta.stats[[tmp.resSubSection]][["TREATMENT"]][["60DPE"]][["Tukey"]] <-
  test.beta.all_metrics$displacement %>%
    dplyr::filter(Treatment  %in% c(
          "A- T- P-",
          # "A- T- P+",
          "A+ T- P-",
          # "A+ T- P+",
          "A- T+ P-",
          # "A- T+ P+",
          "A+ T+ P-"#,
          # "A+ T+ P+"
          )) %>%
    dplyr::filter(Time == 60) %>%
    # Clean up cell value names by removing any strings including and after "__"
    cutCellNames(col = "Beta.Metric", sep = "__") %>%
    # Filter out raw data, leave only normalized data
    dplyr::filter(Data.Type != "Raw") %>%
  # Run Tukey test on GLM models
  run_tukey_glm_beta(., "Displacement", "Beta.Metric", c("Treatment"), 
                family_str = "gaussian",
                group_by_var = NULL) 

Tables

GLM Results - 60 DPE
glm(Displacement ~ Treatment); Prior stressors and no parasite exposure at 60 DPE
term	estimate	std.error	statistic	p.value	p.adj.sig
Bray-Curtis_norm
(Intercept)	−0.060	0.043	−1.398	0.165	ns
A+ T- P-	0.007	0.064	0.109	≥0.25	ns
A- T+ P-	0.043	0.067	0.640	≥0.25	ns
A+ T+ P-	−0.119	0.079	−1.501	0.136	ns
Canberra_norm
(Intercept)	−0.124	0.040	−3.122	0.002	**
A+ T- P-	−0.089	0.059	−1.498	0.137	ns
A- T+ P-	0.099	0.062	1.602	0.112	ns
A+ T+ P-	0.123	0.073	1.680	0.096	ns

ANOVA of GLM - 60 DPE
ANOVA(GLM(Displacement ~ Treatment), type = 2); Prior stressors and no parasite exposure at 60 DPE
term	statistic	df	p.value	sig
Bray-Curtis_norm
Treatment	3.890	3.000	≥0.25	ns
Canberra_norm
Treatment	12.030	3.000	0.007	**

Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); Prior stressors and no parasite exposure at 60 DPE
term	.y.	group1	group2	estimate	std.error	statistic	adj.p.value	Variable
Bray-Curtis_norm
Treatment	Beta.Score	A+ T- P-	A- T- P-	0.007	0.064	0.109	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P-	A- T- P-	0.043	0.067	0.640	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P-	A- T- P-	−0.119	0.079	−1.501	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P-	A+ T- P-	0.036	0.069	0.514	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P-	A+ T- P-	−0.126	0.082	−1.545	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P-	A- T+ P-	−0.162	0.084	−1.931	0.213	Treatment
Canberra_norm
Treatment	Beta.Score	A+ T- P-	A- T- P-	−0.089	0.059	−1.498	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P-	A- T- P-	0.099	0.062	1.602	≥0.25	Treatment
Treatment	Beta.Score	A+ T+ P-	A- T- P-	0.123	0.073	1.680	≥0.25	Treatment
Treatment	Beta.Score	A- T+ P-	A+ T- P-	0.187	0.064	2.921	0.018	Treatment
Treatment	Beta.Score	A+ T+ P-	A+ T- P-	0.212	0.076	2.804	0.025	Treatment
Treatment	Beta.Score	A+ T+ P-	A- T+ P-	0.025	0.078	0.321	≥0.25	Treatment

Significant Results:Pairwise Tukey's HSD, p.adj: Dunnett - 60 DPE
Tukey(Displacement ~ Treatment); All treatments at 60 DPE
term	.y.	group1	group2	estimate	std.error	statistic	adj.p.value	Variable
Canberra_norm
Treatment	Beta.Score	A- T+ P-	A+ T- P-	0.187	0.064	2.921	0.018	Treatment
Treatment	Beta.Score	A+ T+ P-	A+ T- P-	0.212	0.076	2.804	0.025	Treatment

PCoA

Plot

Code

# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressNoParaExp_60DPE"



#### Dispersion --------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]] <-
  run_BetaDispersion(dist.matrix = beta.dist.mat[[tmp.resSubSection]], 
                     var = c("treatment_group"))

##### ANOVA ------------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA"]] <-
  get_HoD_anova(betaDisper = beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]],
                var = c("treatment_group"))

### Table

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA.Table"]] <-
  beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.ANOVA"]] %>%
  
    # Create the table
    dplyr::group_by(Beta.Metric) %>%
    set_GT(var = "p.value", group.by = "Beta.Metric")  %>%
  
    # Title/caption
    gt::tab_header(
      title = "ANOVA: Homogeneity of Dispersion",
      subtitle = "ANOVA(Beta Disperson ~ Treatment); Prior stressors and no parasite exposure at 60 DPE"
    ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = c("p.value", "sig"),
      rows = p.value < 0.05
    )
  )

##### Tukey ------------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey"]] <-
  get_HoD_tukey(betaDisper = beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.model"]],
                var = c("treatment_group")) 

### Table

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey.Table"]] <-
  beta.stats[[tmp.resSubSection]][["TREATMENT"]][["HoD.Tukey"]] %>%
    
    dplyr::mutate(
        group1 = case_when(
        group1 == "Control" ~ "A- T- P-",
        group1 == "Parasite" ~ "A- T- P+",
        group1 == "Antibiotics" ~ "A+ T- P-",
        group1 == "Antibiotics_Parasite" ~ "A+ T- P+",
        group1 == "Temperature" ~ "A- T+ P-",
        group1 == "Temperature_Parasite" ~ "A- T+ P+",
        group1 == "Antibiotics_Temperature" ~ "A+ T+ P-",
        group1 == "Antibiotics_Temperature_Parasite" ~ "A+ T+ P+",
    ),
        group2 = case_when(
        group2 == "Control" ~ "A- T- P-",
        group2 == "Parasite" ~ "A- T- P+",
        group2 == "Antibiotics" ~ "A+ T- P-",
        group2 == "Antibiotics_Parasite" ~ "A+ T- P+",
        group2 == "Temperature" ~ "A- T+ P-",
        group2 == "Temperature_Parasite" ~ "A- T+ P+",
        group2 == "Antibiotics_Temperature" ~ "A+ T+ P-",
        group2 == "Antibiotics_Temperature_Parasite" ~ "A+ T+ P+",
    )) %>%
  
  # Create the table
  dplyr::group_by(Beta.Metric) %>%
  set_GT(var = "adj.p.value", group.by = "Beta.Metric")  %>%
  
  # Title/caption
  gt::tab_header(
    title = "Tukey: Homogeneity of Dispersion",
    subtitle = "Tukey(Beta Disperson ~ Treatment); Prior stressors and no parasite exposure at 60 DPE"
  ) %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = c("adj.p.value", "sig"),
      rows = adj.p.value < 0.05
    )
  )

Tables

ANOVA: Homogeneity of Dispersion
ANOVA(Beta Disperson ~ Treatment); Prior stressors and no parasite exposure at 60 DPE
term	df	sumsq	meansq	statistic	p.value	sig
bray
treatment_group	3.000	0.289	0.096	3.589	0.016	*
Residual	103.000	2.767	0.027	NA	NA	NA
canberra
treatment_group	3.000	0.076	0.025	6.573	<0.001	***
Residual	103.000	0.395	0.004	NA	NA	NA

Tukey: Homogeneity of Dispersion
Tukey(Beta Disperson ~ Treatment); Prior stressors and no parasite exposure at 60 DPE
.y.	term	group1	group2	estimate	conf.low	conf.high	adj.p.value	sig
bray
Distance	treatment_group	A+ T- P-	A- T- P-	−0.036	−0.142	0.070	≥0.25	ns
Distance	treatment_group	A- T+ P-	A- T- P-	−0.064	−0.174	0.046	≥0.25	ns
Distance	treatment_group	A+ T+ P-	A- T- P-	−0.162	−0.294	−0.031	0.009	**
Distance	treatment_group	A- T+ P-	A+ T- P-	−0.028	−0.143	0.087	≥0.25	ns
Distance	treatment_group	A+ T+ P-	A+ T- P-	−0.126	−0.262	0.009	0.077	ns
Distance	treatment_group	A+ T+ P-	A- T+ P-	−0.098	−0.237	0.040	≥0.25	ns
canberra
Distance	treatment_group	A+ T- P-	A- T- P-	−0.056	−0.096	−0.016	0.002	**
Distance	treatment_group	A- T+ P-	A- T- P-	0.010	−0.032	0.051	≥0.25	ns
Distance	treatment_group	A+ T+ P-	A- T- P-	−0.023	−0.072	0.027	≥0.25	ns
Distance	treatment_group	A- T+ P-	A+ T- P-	0.066	0.023	0.109	<0.001	***
Distance	treatment_group	A+ T+ P-	A+ T- P-	0.034	−0.017	0.085	≥0.25	ns
Distance	treatment_group	A+ T+ P-	A- T+ P-	−0.032	−0.085	0.020	≥0.25	ns

CAP

Plot

Code

# Set the subsection for this analysis
tmp.resSubSection <- "PriorStressNoParaExp_60DPE"


#### Capscale ----------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.mod"]] <-
  run_capscale(ps.cleaned_alpha %>% # is the phyloseq object we are looping over
          microViz::ps_filter(Treatment  %in% c(
          "A- T- P-",
          # "A- T- P+",
          "A+ T- P-",
          # "A+ T- P+",
          "A- T+ P-",
          # "A- T+ P+",
          "A+ T+ P-"#,
          # "A+ T+ P+"
          )) %>%
          microViz::ps_filter(Time == 60), 
               dist.matrix = beta.dist.mat[[tmp.resSubSection]], 
               formula_str = "dist ~ Treatment")

##### ADONIS ------------------------------------------------------------------

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS"]] <-
  run_cap_adonis(ps.cleaned_alpha %>% # is the phyloseq object we are looping over
          microViz::ps_filter(Treatment  %in% c(
          "A- T- P-",
          # "A- T- P+",
          "A+ T- P-",
          # "A+ T- P+",
          "A- T+ P-",
          # "A- T+ P+",
          "A+ T+ P-"#,
          # "A+ T+ P+"
          )) %>%
          microViz::ps_filter(Time == 60),
                 dist.matrix = beta.dist.mat[[tmp.resSubSection]], 
                 formula_str = "dist ~ Treatment",
                 by.method = "terms") 

## Warning in tidy.anova(.): The following column names in ANOVA output were not
## recognized or transformed: SumOfSqs, R2

## Warning in tidy.anova(.): The following column names in ANOVA output were not
## recognized or transformed: SumOfSqs, R2

### Table

beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS.Table"]] <-
  beta.stats[[tmp.resSubSection]][["TREATMENT"]][["CAP.ADONIS"]] %>%
    set_GT(var = "p.value", group.by = "Beta.Metric") %>%
  
  # Title/caption
  gt::tab_header(
    title = "ADONIS2",
    subtitle = "adonis2(Beta Distance ~ Treatment); Prior stressors and no parasite exposure at 60 DPE"
  )  %>%
  gt::tab_style(
    style = gt::cell_fill(color = "lightgreen"),
    locations = gt::cells_body(
      columns = c("p.value", "sig"),
      rows = p.value < 0.05
    )
  )

Tables

ADONIS2
adonis2(Beta Distance ~ Treatment); Prior stressors and no parasite exposure at 60 DPE
term	df	SumOfSqs	R2	statistic	p.value	sig
bray
Treatment	3.000	1.688	0.102	3.887	0.001	***
Residual	103.000	14.910	0.898	NA	NA	NA
Total	106.000	16.598	1.000	NA	NA	NA
canberra
Treatment	3.000	3.157	0.086	3.234	0.001	***
Residual	103.000	33.522	0.914	NA	NA	NA
Total	106.000	36.679	1.000	NA	NA	NA